1 Introducción

Un modelo estadístico es una colección de distribuciones de probabilidad indexadas por un parámetro desconocido: $\mathcal{P} = \{p(\mathbf{Y}\mid\theta):\mathbf{Y}\in\mathcal{Y},\theta\in\Theta\}$ donde:

$\theta$ es un parámetro desconocido.
$\Theta$ es el espacio de parámetros, i.e., el conjunto de todos los valores que puede asumir $\theta$ .
$\mathbf{Y}$ es una matriz de adyacencia (sociomatriz).
$\mathcal{Y}$ es el conjunto de todas las posibles matrices de adyacencia.
$p(\mathbf{Y}\mid\theta)$ es la distribución probabilística de $\mathbf{Y}\in\mathcal{Y}$ dado un valor particular de $\theta$ .

Los modelos estadísticos tienen como objetivo:

Investigar el proceso generativo de la red.
Probar la significancia de características predefinidas.
Evaluar las asociaciones entre las variables nodales y relacionales.
Imputar relaciones faltantes.
Predecir relaciones futuras.

A diferencia de los modelos deterministas y algorítmicos, los modelos estadísticos permiten cuantificar la incertidumbre asociada con los procesos.

La riqueza del modelo se deriva en gran medida de cómo especificar $p(\mathbf{Y}\mid\theta)$ .

2 Modelo de grafos aleatorios

El término modelo de grafo aleatorio simple (simple random graph model) o modelo de grafo aleatorio (random graph model) se usa para referirse a un modelo que asume que:

Todos las aristas se forman independientemente unas de las otras.
Todas las aristas tienen una probabilidad común $\theta\in\Theta=(0,1)$ de formarse.

Es decir, las entradas de $\mathbf{Y}$ son independientes e idénticamente distribuidas (iid) de acuerdo con una distribución Bernoulli con parámetro $\theta$ : $y_{i,j}\mid\theta\stackrel{\text{iid}}{\sim} \textsf{Bernoulli}(\theta)$ y por lo tanto $p(\mathbf{Y}\mid\theta) = \prod\theta^{y_{i,j}}(1-\theta)^{1-y_{i,j}} = \theta^{\sum y_{i,j}} (1-\theta)^{\sum(1- y_{i,j})}$ donde la productoria y las sumas se hacen sobre $\{i,j:i\neq j\}$ y $\{i,j:i<j\}$ para redes dirigidas y no dirigidas, respectivamente.

Se utiliza como punto de referencia (baseline) para determinar si un grafo satisface o no alguna propiedad estructural.

Propiedades:

El grado de cada vértice tiene distribución Binomial con parámetros $n-1$ y $\theta$ , donde $n$ es el orden del grafo.
Si $\theta = \frac{k}{n}$ con $k>0$ , entonces la distribución del grado es aproximadamente Poisson con media $k$ cuando $n\rightarrow\infty$ .
Distancias y transitividad bajas.
La distribución de $\mathbf{Y}\mid{\textstyle\sum} y_{i,j} = s,\theta$ es uniforme y no depende de $\theta$ , i.e., el modelo asigna la misma probabilidad a todos los grafos con un número fijo de vértices $n$ y número fijo de aristas $s$ : $p(\mathbf{Y}\mid{\textstyle\sum} y_{i,j} = s,\theta) = \frac{p(\mathbf{Y},{\textstyle\sum y_{i,j}} = s\mid\theta)}{p({\textstyle\sum y_{i,j}} = s\mid\theta)} = \frac{1}{\binom{m}{s}}$ donde $m$ es el número de diadas sobre el cual se hace la sumatoria. Esta distribución se denomina comúnmente como modelo de Erdös–Rényi.

2.1 Ejemplo: Simulación

suppressMessages(suppressWarnings(library(igraph)))
# orden
n <- 100
k <- 2
theta <- k/n
# simulacion (ver tambien sample_gnm)
set.seed(42)
g <- sample_gnp(n = n, p = theta, directed = F, loops = F)
par(mfrow = c(1,2))
igraph_options(vertex.label = NA, vertex.color = 1, vertex.frame.color = "black", vertex.size = 5)
plot(g, layout = layout_in_circle)
plot(g, layout = layout_with_fr)
title(main = "Grafo aleatorio generado con n = 100 y p = 0.02", outer = T, line = -2)

# conectado?
is_connected(g)

## [1] FALSE

# frecuencias de componentes
table(sapply(X = decompose(g), FUN = vcount))

## 
##  1  2  3  4 71 
## 15  2  2  1  1

# distribucion del grado
table(degree(g))/n

## 
##    0    1    2    3    4    5    6 
## 0.15 0.30 0.25 0.17 0.07 0.05 0.01

# grado promedio
mean(degree(g))

## [1] 1.9

# distribucion del grado aprox.
round(dpois(x = 0:6, lambda = k), 3)

## [1] 0.135 0.271 0.271 0.180 0.090 0.036 0.012

# grado promedio aprox.
(n-1)*theta

## [1] 1.98

# distancias promedio
mean_distance(g)

## [1] 5.276511

# clustering coeff.
transitivity(g)

## [1] 0.01639344

2.2 Estimación

¿Cuál valor de $\theta\in\Theta$ que hace $\tilde{\mathbf{Y}}\sim p(\mathbf{Y}\mid\theta)$ sea lo más similar posible a $\mathbf{Y}$ ?

El estimador de máxima verosimilitud (maximum likelihood estimator, MLE) de $\theta$ es el valor $\hat{\theta}\in\Theta$ que maximiza la probabilidad de los datos observados: $p(\mathbf{Y}\mid\hat{\theta})\geq p(\mathbf{Y}\mid\theta)\,\,\text{para todo }\theta\in\Theta\,.$ El MLE de $\theta$ es la densidad de $\mathbf{Y}$ , i.e., $\hat{\theta} = \textsf{den}(\mathbf{Y})$ , dado que la log-verosimilitud está dada por: $\ell(\theta) = \log p(\mathbf{Y}\mid\theta) = \log\theta\textstyle\sum y_{i,j} + \log(1-\theta)\textstyle\sum (1-y_{i,j})\,.$

2.3 Bondad de ajuste

Una prueba de bondad de ajuste interna es una comparación de los datos con el modelo de probabilidad propuesto.

Ingredientes:

Un estadístico de prueba $t$ donde $t:\mathcal{Y}\rightarrow\mathbb{R}$ es una función conocida de los datos.
Una distribución nula $p(t\mid\textsf{M})$ donde $p(t\mid\textsf{M})$ es la distribución de probabilidad de $t$ bajo el modelo estadístico hipotetizado $\textsf{M}$ .
Una medida de discrepancia asociada con la comparación de $t_{\text{obs}} = t(\mathbf{Y})$ con $p(t\mid\textsf{M})$ . Esta medida se denominada valor $p$ y se calcula como $p = \textsf{Pr}(t > t_{\text{obs}}\mid\textsf{M})$ . El valor $p$ se puede aproximar por medio de métodos de Monte Carlo. Valores $p$ que tienden a 0 o a 1 indican una alta discrepancia entre los datos observados y el modelo propuesto.

Una característica se denomina “significativa” si resulta ser inusual o inesperada (valo $p$ extremo) respecto a modelo de grafos aleatorios (e.g., modelo de Erdös–Rényi).

Procedimiento de evaluación del modelo:

Calcular $\hat{\theta}$ a partir del conjunto de datos observado $\mathbf{Y}$ .
Simular $\tilde{\mathbf{Y}}_1,\ldots,\tilde{\mathbf{Y}}_B$ a partir del modelo propuesto $p(\mathbf{Y}\mid\hat{\theta})$ .
Calcular $t_{\text{obs}} = t(\mathbf{Y})$ y $t(\tilde{\mathbf{Y}}_1),\ldots,t(\tilde{\mathbf{Y}}_B)$ .
Comparar $t_{\text{obs}}$ con la distribución empírica de $t(\tilde{\mathbf{Y}}_1),\ldots,t(\tilde{\mathbf{Y}}_B)$ .

Estadísticos de prueba:

Densidad.
Transitividad.
Asortatividad.
Reciprocidad.
Diámetro.
Fracción de vértices aislados, de articulación, o de la componente gigante.
Media, SD, o cualquier característica de una medida de centralidad (grado, closeness, betweenness, eigen).
Media, SD, o cualquier característica de la distancia geodésica.
Cualquier función de los datos observados.

2.4 Ejemplo: Conflictos

Hoff, P. D. (2009). Multiplicative latent factor models for description and prediction of social networks. Computational and mathematical organization theory, 15(4), 261.

Datos de conflictos entre países en los años 90. Y hace referencia a una matriz $\mathbf{Y}=[y_{i,j}]$ en la que $y_{i,j}$ representa el número de conflictos iniciados por el país $i$ hacia el país $j$ .

suppressMessages(suppressWarnings(library(igraph)))
# datos
load("C:/Users/Juan Camilo/Dropbox/UN/networks_2021_2/conflict.RData")
Y <- dat$Y
# simetrizar y binarizar
Y <- 1*(Y+t(Y) > 0)
Y[Y != 0] <- 1
g <- graph_from_adjacency_matrix(adjmatrix = Y, mode = "undirected")
# orden
vcount(g)

## [1] 130

# tamaño
ecount(g)

## [1] 160

# dirigida?
is_directed(g)

## [1] FALSE

# ponderada?
is_weighted(g)

## [1] FALSE

# grafico
igraph_options(vertex.size = 3, edge.color = "darkgray", vertex.frame.color = "black", layout = layout_with_fr)
set.seed(42)
plot(g, vertex.label = NA, main = "Conflictos")

# estimacion de theta MLE
theta_hat <- edge_density(g, loops = FALSE)
theta_hat

## [1] 0.01908169

# log-verosimilitud
n <- vcount(g)
m <- n*(n-1)/2
s <- m*edge_density(g, loops = FALSE)
loglik <- function(theta) s*log(theta) + (m-s)*log(1-theta)
# grafico
par(mfrow = c(1,2))
curve(expr = loglik(x), from = 0, to = 1, lwd = 2, xlab = expression(theta), ylab = "Log-verosimilitud")
abline(v = theta_hat, col = 2, lty = 2)
curve(expr = loglik(x), from = 0, to = 0.04, lwd = 2, xlab = expression(theta), ylab = "Log-verosimilitud")
abline(v = theta_hat, col = 2, lty = 2)

# bondad de ajuste por medio de la transitividad
# no es necesario almacenar los datos simulados, solo los estadisticos
# es posible hacer este computo en paralelo usando la doParallel
t_obs <- transitivity(g)
B <- 1200
t_rep <- NULL
set.seed(42)
for (i in 1:B) {
  g_rep <- sample_gnp(n = n, p = theta_hat, directed = F, loops = F)
  t_rep[i] <- transitivity(g_rep)
}
# grafico
hist(x = t_rep, freq = F, col = "gray90", border = "gray90", xlim = c(0,0.2), xlab = "Transitividad", ylab = "Densidad", main = "")
abline(v = t_obs, col = 4, lty = 2)

# valor p
mean(t_rep > t_obs)

## [1] 0

2.5 Ejemplo: Áreas cerebrales visotáctiles

Negyessy L., Nepusz T., Kocsis L., Bazso F.: Prediction of the main cortical areas and connections involved in the tactile function of the visual cortex by network analysis. European Journal of Neuroscience, 23(7): 1919-1930, 2006.

Grafo de las áreas cerebrales visuotáctiles del mono macaco. El modelo consta de 45 áreas y 463 conexiones dirigidas.

suppressMessages(suppressWarnings(library(igraphdata)))
data(macaque)
# orden
vcount(macaque)

## [1] 45

# tamaño
ecount(macaque)

## [1] 463

# dirigida?
is_directed(macaque)

## [1] TRUE

# ponderada?
is_weighted(macaque)

## [1] FALSE

# grafico
igraph_options(vertex.size = 6, edge.arrow.size = 0.3, edge.color = "darkgray", vertex.frame.color = "black", layout = layout_with_fr)
set.seed(42)
plot(macaque, vertex.label = NA, main = "Áreas cerebrales visotáctiles")

# datos
data(macaque)
nv <- vcount(macaque)
ne <- ecount(macaque)
# numero de replicas
ntrials <- 1000
# para almacenar numero de comunidades
cl.rg <- NULL
apl.rg <- NULL
# simulacon
set.seed(42)
for (i in 1:ntrials) {
   g.rg <- sample_gnm(nv, ne, directed=TRUE)
   cl.rg[i] <- transitivity(g.rg)
   apl.rg[i] <- mean_distance(g.rg)
}

# grafico
par(mfrow = c(1,2))
# transitividad
hist(x = cl.rg, freq = F, col = "gray90", border = "gray90", xlim = c(0.35,0.55), xlab = "Transitividad", ylab = "Densidad", main = "")
abline(v = transitivity(macaque), col = 4, lty = 2)
# distancia promedio
hist(x = apl.rg, freq = F, col = "gray90", border = "gray90", xlim = c(1.8,2.2), xlab = "Distancia prom.", ylab = "Densidad", main = "")
abline(v = mean_distance(macaque), col = 4, lty = 2)

3 Modelo de grafos aleatorio generalizado

El modelo de Erdös–Rényi se puede generalizar de manera sencilla:

Definir la colección de todos los grafos con un orden fijo $n$ que posean una característica determinada.
Asignar la misma probabilidad a cada uno de los grafos de esta colección.

La característica más común es la de una secuencia de grado fija $\{d_{(1)},\ldots,d_{(n)}\}$ (en forma ordenada).

Los grafos de la colección con un número fijo de vértices y una secuencia de grado fija tienen el mismo número de aristas $s$ , debido a que el grado medio es $\bar{d} = 2s/n$ . Por lo tanto, esta colección está estrictamente contenida dentro de la colección de grafos generados a partir del modelo de Erdös–Rényi.

Por otro lado, es importante tener en cuenta que todas las demás características pueden variar libremente en la medida que lo permita la secuencia de grados elegida.

Las rutinas tipo sample_[...] en igraph ofrecen más alternativas para generar modelos de grafos aleatorios generalizados.

3.1 Ejemplo: Simulación

# secuencia de grado
degs <- c(2,2,2,2,3,3,3,3)
# simulacion
set.seed(42)
g1 <- sample_degseq(degs, method = "vl")
g2 <- sample_degseq(degs, method = "vl")
# grafico
par(mfrow = c(1,2))
igraph_options(vertex.size = 8, edge.color = "darkgray", vertex.frame.color = "black", layout = layout_with_fr)
plot(g1, vertex.label=NA)
plot(g2, vertex.label=NA)

# se satisface la condicion?
all(degree(g1) == degs)

## [1] TRUE

all(degree(g2) == degs)

## [1] TRUE

# isomorfos?
isomorphic(g1, g2)

## [1] FALSE

# numero de aristas
c(ecount(g1), ecount(g2))

## [1] 10 10

# grado medio
c(mean(degs), 2*ecount(g1)/vcount(g1), 2*ecount(g2)/vcount(g2))

## [1] 2.5 2.5 2.5

3.2 Ejemplo: Simulación

# Power-law degree distribution
# corregir si la suma de la secuencia de grado es impar 
set.seed(42)
degs <- sample(x = 1:100, size = 100, replace = TRUE, prob = (1:100)^-2)
# Exponential degree distribution
# degs <- sample(1:100, 100, replace=TRUE, prob=exp(-0.5*(1:100)))
if (sum(degs) %% 2 != 0) { degs[1] <- degs[1] + 1 }
g <- sample_degseq(degs, method = "vl")
# se satisface la condicion?
all(degree(g) == degs)

## [1] TRUE

# grafico
igraph_options(vertex.size = 6, edge.color = "darkgray", vertex.frame.color = "black", layout = layout_with_fr)
plot(g, vertex.label=NA)

# grado
d <- degree(g)
# distribucion de grado
dd <- degree_distribution(g)
# grafico
par(mfrow = c(1,2))
plot(NA, NA, type = "n", xlim = c(0,50), ylim = c(0,0.17), xlab = "Grado", ylab = "Densidad", main = "Distribución de grado")
hist(d, freq = F, col = "gray90", border = "gray50", add = T)
ind <- dd != 0
plot((0:max(d))[ind], dd[ind], log = "xy", col = "blue", xlab = "Log-grado", ylab = "Log-densidad", main = "Distribución de grado (log-log)")

3.3 Ejemplo: Yeast

Von Mering, C., Krause, R., Snel, B., Cornell, M., Oliver, S. G., Fields, S., & Bork, P. (2002). Comparative assessment of large-scale data sets of protein–protein interactions. Nature, 417(6887), 399-403.

Interacciones entre proteínas que tienen una confianza “alta” y “media”. Las interacciones proteína-proteína prometen revelar aspectos de la compleja red reguladora que subyace a la función celular.

# datos
suppressMessages(suppressWarnings(library(igraphdata)))
data(yeast)
# distribucion de grado real
degs <- degree(yeast)
# simulacion
fake.yeast <- sample_degseq(degs, method = c("vl"))
# se satisface la condicion?
all(degree(yeast) == degree(fake.yeast))

## [1] TRUE

# se mantienen otras propiedades de la red?
# diametro
diameter(yeast)

## [1] 15

diameter(fake.yeast)

## [1] 8

# transitividad
 transitivity(yeast)

## [1] 0.4686178

transitivity(fake.yeast)

## [1] 0.03809097

3.4 Ejemplo: Zachary

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

Los nodos representan a los miembros de un club de karate observado durante aproximadamente 2 años y los enlaces que conectan dos nodos indican interacciones sociales entre los dos miembros correspondientes.

# datos 
data(karate)
nv <- vcount(karate)
ne <- ecount(karate)
degs <- degree(karate)
# numero de replicas
ntrials <- 1000
# para almacenar numero de comunidades
num.comm.rg <- NULL 
num.comm.grg <- NULL
# simulacion
set.seed(42)
for(i in 1:ntrials){
  # modelo de grafos aleatorio simple
  g.rg <- sample_gnm(n = nv, m = ne, directed = F, loops = F)
  c.rg <- cluster_fast_greedy(graph = g.rg)
  num.comm.rg[i] <- length(c.rg)
  # modelo de grafos aleatorio generalizado
  g.grg <- sample_degseq(degs, method = "vl")
  c.grg <- cluster_fast_greedy(graph = g.grg)
  num.comm.grg[i] <- length(c.grg)
}

# distribucion empirica
rslts <- c(num.comm.rg, num.comm.grg)
indx <- c(rep(0, ntrials), rep(1, ntrials))
counts <- table(indx, rslts)/ntrials
barplot(counts, beside = T, col = c("blue","red"), xlab = "Número de comunidades", ylab = "Frecuencia relativa", legend = c("Modelo SRG","Modelo GSRG"))

# valor observado
length(cluster_fast_greedy(graph = karate))

## [1] 3

Es muy probable que existan mecanismos adicionales en la generación de la dináica social de la red de karate, más allá de la simple densidad y distribución de las interacciones sociales en la red.

4 Modelo de mundo pequeño

Muchas redes muestran altos niveles de agrupación, pero pequeñas distancias entre los nodos. Tal comportamiento no se puede reproducir mediante un modelo de grafos aleatorios (pueden reproducir pequeñas distancias pero con bajos niveles de agrupación).

En el modelo de mundo pequeño (small-world model) los $n$ vértices se disponen de manera periódica en un arreglo (lattice) de $d$ dimensiones, y se enlazan a sus $r$ vecinos en cada dirección del arreglo. Luego, un extremo de cada arista, independientemente y con probabilidad $p$ , se mueve para incidir en otro vértice elegido de manera uniforme.

Volver a “cablear” aleatoriamente una fracción pequeña de aristas reduce notablemente la distancia entre los vértices, manteniendo un nivel de agrupamiento igualmente alto.

4.1 Ejemplo: Simulación

# simulacion
# d = 1, n = 12, r = 2
# varios valores de p
par(mfrow = c(3,3))
set.seed(42)
for (p in seq(from = 0, to = 1, len = 9)) {
  g <- sample_smallworld(dim = 1, size = 12, nei = 2, p = p, loops = F, multiple = F)
  plot(g, layout = layout_in_circle, vertex.label = NA, vertex.size = 12, edge.color = "darkgray", main = paste0("p = ", round(p, 2)))
}

# simulacion
# d = 1, n = 100, r = 5, p = 0
set.seed(42)
g1 <- sample_smallworld(dim = 1, size = 100, nei = 5, p = 0, loops = F, multiple = F)
set.seed(42)
g2 <- sample_smallworld(dim = 1, size = 100, nei = 5, p = 0.05, loops = F, multiple = F)
par(mfrow = c(1,2))
plot(g1, layout = layout_in_circle, vertex.label = NA, vertex.size = 4, edge.color = "darkgray", main = "p = 0")
plot(g2, layout = layout_in_circle, vertex.label = NA, vertex.size = 4, edge.color = "darkgray", main = "p = 0.05")

# diametro
c(diameter(g1), diameter(g2))

## [1] 10  5

# distancia promedio
c(mean_distance(g1), mean_distance(g2))

## [1] 5.454545 2.689697

# transitividad
c(transitivity(g1), transitivity(g2))

## [1] 0.6666667 0.4953704

# simulacion
steps <- seq(from = -4, to = -0.5, by = 0.1)
len <- length(steps)
cl  <- NULL
apl <- NULL
ntrials <- 100
set.seed(42)
for (i in 1:len) {
 cltemp <- numeric(ntrials)
 apltemp <- numeric(ntrials)
 for (j in (1:ntrials)) {
   g <- sample_smallworld(dim = 1, size = 1000, nei = 10, p = 10^steps[i], loops = F, multiple = F)
   cltemp[j] <- transitivity(g)
   apltemp[j] <- mean_distance(g)
 }
 cl[i]  <- mean(cltemp)
 apl[i] <- mean(apltemp)
}
# grafico
plot(steps, cl/max(cl), ylim = c(0, 1), lwd = 3, type = "l", col = "blue", xlab = expression(log[10](p)), ylab = "Coeficientes")
lines(steps, apl/max(apl), lwd = 3, col = "red")
legend("bottomleft", legend = c("Transitividad","Distancia media"), col = c("blue","red"), lwd = 3, bty = "n")

5 Modelo de fijación preferencial

Muchas redes crecen o evolucionan con el tiempo.

World Wide Web.
Redes de citas científicas.
Redes biológicas.

Hay varios mecanismos para especificar como una red evoluciona a través del tiempo. Uno de tales mecanismos es el modelo de fijación preferencial (preferential attachment model) o modelo de Barabasi-Albert, diseñado para considerar tanto el crecimiento como el principio de que “los ricos se hacen más ricos” (‘the rich get richer’). Reproducir distribuciones de grado donde los nodos populares tienden a acumular un número cada vez mayor de enlaces a medida que pasa el tiempo.

Estado inicial: grafo $G^{(0)}$ con $n^{(0)}$ vértices y $s^{(0)}$ aristas.
En el estado $t=1,2,\ldots$ se genera un nuevo grafo $G^{(t)}$ a partir de $G^{(t-1)}$ agregando un nuevo vértice de grafo $m\geq 1$ de tal forma que las nuevas $m$ aristas se fijen a $m$ vértices diferentes en $G^{(t-1)}$ con probabilidad proporcional al grado, $p_i\propto \left(d_i^{(t-1)}\right)^\beta + \gamma$ , para $i = 1,\ldots,n^{(t-1)}$ , con $\beta,\gamma >0$ (comúnmente $\beta = \gamma = 1$ ).

En cada etapa, $m$ vértices existentes están conectados a un nuevo vértice de una manera preferencial a aquellos con grados más altos. Después de $t$ iteraciones $G^{(t)}$ tendrá $n^{(t)} = n^{(0)} + t$ vértices y $s^{(t)} = s^{(0)} + tm$ aristas. Los vértices de grado alto deberían emerger gradualmente a medida que $t$ aumenta.