Processing math: 36%

1 Introducción

Un modelo estadístico se define como una colección de distribuciones de probabilidad parametrizadas por un parámetro de dimensión finita desconocido:
P={p(Yθ):YY,θΘ}
donde:

El propósito principal de los modelos estadísticos es:

Además, estos modelos permiten cuantificar la incertidumbre inherente a dichos procesos.

La especificación de p(Yθ) es clave para determinar la riqueza del modelo, lo cual puede abordarse mediante diversas aproximaciones, entre ellas:

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 refiere a un modelo en el que:

En este modelo, las entradas de la matriz de adyacencia Y son independientes e idénticamente distribuidas (iid) según una distribución Bernoulli con parámetro θ:
yi,jθiidBernoulli(θ),
lo que implica que la función de probabilidad conjunta es:
p(Yθ)=θyi,j(1θ)1yi,j=θyi,j(1θ)(1yi,j),
donde el producto y las sumas se toman sobre {i,j:ij} para redes dirigidas y {i,j:i<j} para redes no dirigidas.

Este modelo sirve como un punto de referencia (baseline) para analizar si un grafo satisface ciertas propiedades locales o estructurales.

2.0.1 Propiedades del Modelo

  1. Distribución de grado:
    • El grado de cada vértice sigue una distribución Binomial con parámetros n1 y θ, donde n=|V| es el número de nodos del grafo.
    • Si θ=kn con k>0, la distribución del grado converge aproximadamente a una distribución Poisson con media k cuando n.
  2. Distancias y transitividad:
    • El modelo típicamente produce grafos con bajas distancias promedio y baja transitividad.
  3. Distribución condicional uniforme:
    • Para un número fijo de aristas s, la distribución de Yyi,j=s,θ es uniforme y no depende de θ. Esto implica que el modelo asigna la misma probabilidad a todos los grafos con n nodos y s aristas:
      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{q}{s}},
      donde q es el número de díadas en las que se define la sumatoria.

Este modelo, conocido comúnmente como modelo de Erdős–Rényi, es fundamental para entender las propiedades probabilísticas de los grafos aleatorios.

2.1 Ejemplo: Simulación

suppressMessages(suppressWarnings(library(igraph)))
# orden
n <- 100
k <- 2
theta <- k/n
# simulación (ver también sample_gnm)
set.seed(42)
g <- sample_gnp(n = n, p = theta, directed = F, loops = F)
igraph_options(vertex.label = NA, edge.color = "gray40", vertex.color = 1, vertex.frame.color = 1, vertex.size = 6)
set.seed(42)
plot(g, layout = layout_with_fr, main = "Grafo aleatorio generado con n = 100 y p = 0.02")

# 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
# distribución 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
# distribución 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
mean(degree(g))
## [1] 1.9
# grado promedio aprox.
(n-1)*theta
## [1] 1.98
# distancia promedio
mean_distance(g)
## [1] 5.276511
# coeficiente de agrupamiento
transitivity(g)
## [1] 0.01639344

2.2 Ejemplo: Simulación (cont.)

suppressMessages(suppressWarnings(library(igraph)))
# orden
n <- 100
k <- 2
theta <- k/n
# simulación
B  <- 1000
dg <- NULL
dp <- NULL
ca <- NULL
set.seed(42)
for (b in 1:B) {
  # grafo
  g  <- sample_gnp(n = n, p = theta, directed = F, loops = F)
  # estadísticos
  dg <- rbind(dg, table(factor(x = degree(g), levels = 0:9))/n)
  dp[b] <- mean_distance(g)
  ca[b] <- transitivity(g)
}
# distribución del grado
par(mar = c(3,3,1,1), mgp = c(1.75,0.75,0))
dg_q <- apply(X = dg, MARGIN = 2, FUN = quantile, probs = c(0.025,0.5,0.975))
plot(NA, NA, xlim = c(0,9), ylim = range(dg_q), xlab = "Grado", ylab = "Densidad")
abline(v = 0:9, col = "lightgray")
segments(x0 = 0:9, y0 = dg_q[1,], x1 = 0:9, y1 = dg_q[3,])
lines(x = 0:9, y = dg_q[2,], type = "p", pch = 16)
lines(x = 0:9, y = dpois(x = 0:9, lambda = k), col = 4, type = "p", pch = 17)
legend("topright", legend = c("Sim.", "Pois."), col = c(1,4), pch = c(16,17))

par(mfrow = c(1,2), mar = c(3,3,1,1), mgp = c(1.75,0.75,0))
hist(x = dp, freq = F, col = "gray90", border = "gray90", xlab = "Distancia prom.", ylab = "Densidad", main = "Distancia prom.")
abline(v = quantile(x = dp, probs = c(0.025,0.5,0.975)), col = c(2,4,2), lty = c(2,4,2) ,lwd = c(2,1,2))
hist(x = ca, freq = F, col = "gray90", border = "gray90", xlab = "Transitividad", ylab = "Densidad", main = "Transitividad")
abline(v = quantile(x = ca, probs = c(0.025,0.5,0.975)), col = c(2,4,2), lty = c(2,4,2) ,lwd = c(2,1,2))
legend("topright", legend = c("P 50", "IC 95"), col = c(2,4), lty = c(2,4), lwd = 2)

2.3 Estimación

¿Cuál es el valor de \theta \in \Theta que hace que \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.

Formalmente:
p(\mathbf{Y} \mid \hat{\theta}) \geq p(\mathbf{Y} \mid \theta) \quad \text{para todo } \theta \in \Theta.

En este caso, el MLE de \theta corresponde a la densidad de la matriz \mathbf{Y}, es decir,
\hat{\theta} = \textsf{den}(\mathbf{Y}),
donde \textsf{den}(\mathbf{Y}) es la proporción de entradas activas (aristas presentes) en \mathbf{Y}.

La log-verosimilitud, que sirve como función objetivo para encontrar \hat{\theta}, 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}),
donde las sumas recorren las entradas relevantes de \mathbf{Y}.

2.4 Bondad de ajuste

Una prueba de bondad de ajuste interna evalúa la concordancia entre los datos observados y el modelo probabilístico propuesto.

2.5 Elementos clave

  • Estadístico de prueba:
    Una función real t: \mathcal{Y} \to \mathbb{R} que resume características relevantes de los datos observados.

  • Distribución nula:
    La distribución p(t \mid \textsf{M}), que describe el comportamiento esperado de t bajo el modelo estadístico \textsf{M}.

  • Medida de discrepancia:
    Se utiliza el valor p, que mide la probabilidad de observar valores del estadístico de prueba mayores que el valor observado t_{\text{obs}} = t(\mathbf{Y}) bajo el modelo \textsf{M}: p = \textsf{Pr}(t > t_{\text{obs}} \mid \textsf{M}). Este valor p puede aproximarse utilizando métodos de Monte Carlo. Valores extremos (cercanos a 0 o 1) indican una discrepancia significativa entre los datos y el modelo propuesto.

Una característica se considera significativa si su valor p es extremo, lo que sugiere que el modelo no explica bien dicha característica (por ejemplo, en modelos de grafos aleatorios como el de Erdős–Rényi).

2.6 Procedimiento para Evaluar el Modelo

  1. Estimar parámetros:
    Calcular \hat{\theta} a partir de los datos observados \mathbf{Y}.

  2. Simular datos:
    Generar B muestras simuladas \tilde{\mathbf{Y}}_1, \ldots, \tilde{\mathbf{Y}}_B a partir del modelo p(\mathbf{Y} \mid \hat{\theta}).

  3. Calcular el estadístico:
    Evaluar el estadístico de prueba para los datos observados t_{\text{obs}} = t(\mathbf{Y}) y para las simulaciones t(\tilde{\mathbf{Y}}_1), \ldots, t(\tilde{\mathbf{Y}}_B).

  4. Comparar resultados:
    Comparar t_{\text{obs}} con la distribución empírica de t(\tilde{\mathbf{Y}}_1), \ldots, t(\tilde{\mathbf{Y}}_B) para calcular el valor p.

2.7 Ejemplos de Estadísticos de Prueba

  • Propiedades estructurales:
    • Densidad.
    • Transitividad.
    • Asortatividad.
    • Reciprocidad.
    • Diámetro.
  • Características nodales y de la red:
    • Fracción de vértices aislados, de articulación o en la componente gigante.
    • Media, desviación estándar (SD) u otras estadísticas relacionadas con medidas de centralidad.
    • Media, SD o características de las distancias geodésicas.
  • Funciones generales:
    Cualquier función de los datos observados que capture propiedades relevantes.

2.8 Ejemplo: Conflictos

Los datos hacen corresponden a los conflictos entre países durante la década de 1990. La matriz \mathbf{Y} = [y_{i,j}] representa las relaciones de conflicto entre pares de países, donde cada entrada y_{i,j} indica el número de conflictos iniciados por el país i hacia el país j.

Hoff, P. D. (2009). Multiplicative latent factor models for description and prediction of social networks. Computational and Mathematical Organization Theory, 15(4), 261–272.

suppressMessages(suppressWarnings(library(igraph)))
# datos
load("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
# Cargar las bibliotecas necesarias
suppressMessages(suppressWarnings(library(igraph)))
suppressMessages(suppressWarnings(library(ggraph)))
suppressMessages(suppressWarnings(library(tidygraph)))
# Crear un grafo no dirigido
g <- graph_from_adjacency_matrix(adjmatrix = Y, mode = "undirected")
# Convertir a un formato compatible con ggraph/tidygraph
g_tidy <- as_tbl_graph(g)
# Calcular el grado de los nodos
g_tidy <- g_tidy %>%
  mutate(degree = centrality_degree()) 
# Visualizar la red con ggraph
set.seed(123)  
ggraph(g_tidy, layout = "fr") +
  geom_edge_link(aes(alpha = 0.5), color = "gray", width = 0.3) +
  geom_node_point(aes(size = degree), color = "skyblue") +
  theme_void() +
  labs(title = "Red de conflictos entre países (1990s)") +
  theme(legend.position = "none")

# estimación de theta MLE
theta_hat <- edge_density(g, loops = FALSE)
theta_hat
## [1] 0.01908169
mean(Y[lower.tri(Y)])
## [1] 0.01908169
mean(Y[upper.tri(Y)])
## [1] 0.01908169
n <- dim(Y)[1]
sum(Y)/(n*(n-1))
## [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)
# gráfico
par(mfrow = c(1,2), mar = c(3,3,1,1), mgp = c(1.75,0.75,0))
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 <- 1000
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)
}
# gráfico
hist(x = t_rep, freq = F, col = "gray90", border = "gray90", xlim = c(0,0.2), xlab = "Transitividad", ylab = "Densidad", main = "Conflictos")
abline(v = t_obs, col = 4, lty = 2)

# valor p
mean(t_rep > t_obs)
## [1] 0

3 Validación Cruzada

La validación cruzada consiste en evaluar un modelo particionando el conjunto de díadas D en grupos para medir su capacidad predictiva.

El procedimiento es el siguiente:

  1. Dividir el conjunto de díadas D en L grupos (folds) C_1, \ldots, C_L, procurando que los grupos sean de tamaño similar.
  2. Para cada grupo \ell = 1, \ldots, L:
    • Definir díadas de entrenamiento: Utilizar como datos de entrenamiento las díadas en D - C_\ell, es decir, todas las díadas excepto las del grupo de prueba actual.
    • Definir díadas de prueba: Usar las díadas en C_\ell como datos de prueba.
    • Ajustar el modelo: Entrenar el modelo utilizando únicamente las díadas de entrenamiento.
    • Calcular probabilidades de interacción: Generar las probabilidades de interacción para las díadas de prueba con el modelo ajustado.
    • Evaluar el modelo: Comparar las probabilidades de interacción predichas con los valores observados de las díadas de prueba.

3.1 Matriz de confusión

La matriz de confusión es una herramienta clave para evaluar modelos de clasificación, comparando las categorías reales (positivas o negativas) con las predicciones del modelo. Los valores principales son:

  • Verdaderos Positivos (TP): Casos positivos correctamente identificados por el modelo.
  • Verdaderos Negativos (TN): Casos negativos correctamente identificados por el modelo.
  • Falsos Positivos (FP): Casos negativos que el modelo predijo como positivos (errores tipo I).
  • Falsos Negativos (FN): Casos positivos que el modelo predijo como negativos (errores tipo II).

A partir de estos valores, se calculan métricas clave como:

  • Sensibilidad (recall): mide la capacidad de identificar correctamente los positivos: \text{Sensibilidad} = \frac{\text{TP}}{\text{TP} + \text{FN}}

  • Especificidad: mide la capacidad de identificar correctamente los negativos: \text{Especificidad} = \frac{\text{TN}}{\text{TN} + \text{FP}}

  • Precisión: mide la proporción de predicciones positivas que son correctas: \text{Precisión} = \frac{\text{TP}}{\text{TP} + \text{FP}}

  • Exactitud: mide la proporción total de predicciones correctas: \text{Exactitud} = \frac{\text{TP} + \text{TN}}{\text{TP} + \text{TN} + \text{FP} + \text{FN}}

Estas métricas ofrecen una evaluación integral del desempeño del modelo, permitiendo identificar fortalezas y áreas de mejora.

3.2 Curva ROC

Una curva ROC (Receiver Operating Characteristic curve) es una representación gráfica que muestra la relación entre la Tasa de Verdaderos Positivos (Sensibilidad) y la Tasa de Falsos Positivos (1 - Especificidad) para un clasificador binario, al modificar el umbral de decisión.

Esta curva permite evaluar el desempeño del modelo a través de diferentes niveles de tolerancia a errores, ayudando a identificar el balance óptimo entre sensibilidad y especificidad.

3.3 Ejemplo: Conflictos (cont.)

suppressMessages(suppressWarnings(library(igraph)))
# datos
load("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")

Conformación de folds

# orden
(n <- vcount(g))
## [1] 130
# numero de diadas
(M <- n*(n-1)/2)
## [1] 8385
# numero de folds
(L <- 10)
## [1] 10
# conformación de folds
set.seed(42)
fold_index_vec <- sample(x = 1:L, size = M, replace = T)
fold_index_mat <- matrix(data = 0, nrow = n, ncol = n)
fold_index_mat[lower.tri(fold_index_mat)] <- fold_index_vec
# distribución de folds
(tab <- table(fold_index_vec))
## fold_index_vec
##   1   2   3   4   5   6   7   8   9  10 
## 854 831 815 838 849 800 816 800 846 936
# viz de folds a traves de la matriz de adyacencia
corrplot::corrplot(corr = fold_index_mat/L, col.lim = c(0,1), method = "color", tl.pos = "n", cl.pos = "n")

# distribución de enlaces por fold
y <- Y[lower.tri(Y)]
tmp <- NULL
for (l in 1:L)
  tmp <- rbind(tmp, table(y[fold_index_vec == l])/tab[l])
round(tmp, 3)
##           0     1
##  [1,] 0.980 0.020
##  [2,] 0.986 0.014
##  [3,] 0.984 0.016
##  [4,] 0.969 0.031
##  [5,] 0.987 0.013
##  [6,] 0.984 0.016
##  [7,] 0.977 0.023
##  [8,] 0.978 0.022
##  [9,] 0.979 0.021
## [10,] 0.986 0.014

Probabilidades de interacción en cada fold

# validación cruzada
IP <- vector(mode = "list", L)
B <- 1000
set.seed(123)
for (l in 1:L) {
  # datos de entrenamiento
  y_train <- y  
  y_train[fold_index_vec == l] <- NA
  # ajuste del modelo
  theta_hat <- mean(y_train, na.rm = T)
  # predecir
  n_test <- tab[l]
  inter_prob <- rep(0, n_test)
  for (b in 1:B)
    inter_prob <- inter_prob + rbinom(n = n_test, size = 1, prob = theta_hat)/B
  IP[[l]] <- inter_prob
}

Curvas ROC

# curvas ROC y AUCs
aucs <- NULL
plot(NA, NA, xlim = c(0,1), ylim = c(0,1), xlab = "Tasa Falsos Positivos", ylab = "Tasa verdaderos positivos", main = "Curva ROC")
abline(a = 0, b = 1, col = "gray", lwd = 2)
for (l in 1:L) {
  # datos de prueba
  y_test <- y[fold_index_vec == l]
  # rendimiento
  pred <- ROCR::prediction(predictions = IP[[l]], labels = y_test)
  perf <- ROCR::performance(prediction.obj = pred, measure = "tpr", x.measure = "fpr")
  # ROC
  lines(x = perf@x.values[[1]], y = perf@y.values[[1]], type = "l", col = 2)
  # AUC
  perf <- ROCR::performance(prediction.obj = pred, measure = "auc")
  aucs[l] <- perf@y.values[[1]]
}

Áreas bajo la curva ROC

El AUC (Area Under the Curve) en una curva ROC mide el desempeño general de un clasificador binario como el área bajo la curva que relaciona la Tasa de Verdaderos Positivos (Sensibilidad) y la Tasa de Falsos Positivos (1 - Especificidad) para diferentes umbrales de decisión.

  • Rango: Entre 0 y 1.
    • AUC = 1: Modelo perfecto.
    • AUC = 0.5: Modelo sin capacidad discriminativa (equivalente a azar).

El AUC permite comparar clasificadores y evaluar su capacidad para distinguir entre clases, incluso en datos desbalanceados. Cuanto más cercano a 1, mejor el rendimiento.

# AUCs
round(aucs, 4)
##  [1] 0.3807 0.4939 0.5775 0.4189 0.4876 0.5928 0.3973 0.4791 0.5258 0.6258
# AUC promedio
round(mean(aucs), 4)
## [1] 0.4979
# AUC CV
round(sd(aucs)/mean(aucs), 2)
## [1] 0.17

4 Referencias