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∣θ):Y∈Y,θ∈Θ}
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:
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∣θiid∼Bernoulli(θ),
lo que implica que la función de probabilidad conjunta es:
p(Y∣θ)=∏θyi,j(1−θ)1−yi,j=θ∑yi,j(1−θ)∑(1−yi,j),
donde el producto y las sumas se toman sobre {i,j:i≠j} 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.
Este modelo, conocido comúnmente como modelo de Erdős–Rényi, es fundamental para entender las propiedades probabilísticas de los grafos aleatorios.
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
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)
¿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}.
Una prueba de bondad de ajuste interna evalúa la concordancia entre los datos observados y el modelo probabilístico propuesto.
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).
Estimar parámetros:
Calcular \hat{\theta} a partir de
los datos observados \mathbf{Y}.
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}).
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).
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.
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
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:
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:
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.
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.
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
(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
# 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 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]]
}
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.
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