Modelo de Sociabilidad

Modelo

Considere una red binaria no dirigida representada mediante la matriz de adyacencia \(\mathbf{Y} = [y_{i,j}]\), donde cada elemento \(y_{i,j} \in \{0,1\}\) indica la presencia (\(y_{i,j} = 1\)) o ausencia (\(y_{i,j} = 0\)) de una conexión entre los \(n\) individuos que conforman la red.

El modelo para la probabilidad de conexión entre dos nodos se define como:
\[ y_{i,j} \mid \theta_{i,j} \overset{\text{ind}}{\sim} \textsf{Ber}(\theta_{i,j}), \quad \text{para } i < j, \]
con:
\[ \theta_{i,j} = \Phi(\mu + \delta_i + \delta_j), \]
donde \(\Phi(\cdot)\) denota la función de distribución acumulada de una Normal estándar. Esta función se utilizada como función de enlace para garantizar que las probabilidades \(\theta_{i,j,k}\) permanezcan dentro del intervalo \((0, 1)\).

Para los parámetros, se asumen las siguientes distribuciones previas:
\[ \mu \sim \textsf{N}(0, \sigma^2), \quad \delta_i \overset{\text{iid}}{\sim} \textsf{N}(0, \tau^2), \quad \sigma^2 \sim \textsf{GI}(a_\sigma, b_\sigma), \quad \tau^2 \sim \textsf{GI}(a_\tau, b_\tau). \]

El conjunto de parámetros del modelo es:
\[ \Theta = \{\mu, \delta_1, \delta_2, \dots, \delta_n, \sigma^2, \tau^2\}, \] donde:

• \(\mu\): Representa el efecto promedio global de conectividad en la red; valores altos indican una mayor tendencia general a formar conexiones.
  
• \(\delta_i\) (\(i = 1, \dots, n\)): Captura la sociabilidad del nodo \(i\); valores positivos indican una mayor propensión individual a conectarse, mientras que valores negativos reflejan una menor tendencia.
  
• \(\sigma^2\): Varianza de \(\mu\), mide la incertidumbre en el efecto global de conectividad.
  
• \(\tau^2\): Varianza de \(\delta_i\), refleja la heterogeneidad en las tendencias de conexión entre nodos.

Los \(\delta_i\) capturan diferencias individuales en conectividad, permitiendo modelar patrones de centralidad y aislamiento en la red.

En total, el modelo tiene \(n + 3\) parámetros: \(n\) efectos de sociabilidad, más \(\mu\), \(\sigma^2\), y \(\tau^2\).

Los hiperparámetros del modelo son:
\[ \{a_\sigma, b_\sigma, a_\tau, b_\tau\}, \]
donde \(a_\sigma, b_\sigma\) controlan la varianza de \(\mu\), y \(a_\tau, b_\tau\) regulan la variabilidad de \(\delta_i\).

Estimación

En el modelo propuesto, las distribuciones condicionales completas no tienen forma cerrada debido a que la función de enlace introduce una relación no lineal entre los parámetros y las observaciones.

Uso de variables latentes auxiliares

Para abordar esta limitación, se puede introducir variables auxiliares latentes \(z_{i,j}\), tales que:
\[ y_{i,j} = \begin{cases} 1 & \text{si } z_{i,j} > 0, \\ 0 & \text{si } z_{i,j} \leq 0, \end{cases} \] donde: \[ z_{i,j} \mid \mu, \delta_i, \delta_j \sim \textsf{N}(\mu + \delta_i + \delta_j, 1). \]

La introducción de las variables latentes \(z_{i,j}\) mantiene el modelo original porque estas variables son consistentes con la probabilidad subyacente del modelo inicial. Al integrar \(z_{i,j}\) sobre su distribución condicional, recuperamos la probabilidad original de \(y_{i,j}\):

1. Dado que \(z_{i,j} \mid \mu, \delta_i, \delta_j \sim \textsf{N}(\mu + \delta_i + \delta_j, 1)\), la probabilidad de que \(y_{i,j} = 1\) es equivalente a la probabilidad de que \(z_{i,j} > 0\):
   \[ \textsf{P}(y_{i,j} = 1 \mid \mu, \delta_i, \delta_j) = \textsf{P}(z_{i,j} > 0) = \int_0^\infty \textsf{N}(z_{i,j} \mid \mu + \delta_i + \delta_j, 1) \, \textsf{d}z_{i,j}. \]

2. Evaluando esta integral, obtenemos: \[ \textsf{P}(y_{i,j} = 1 \mid \mu, \delta_i, \delta_j) = \Phi(\mu + \delta_i + \delta_j), \] donde \(\Phi(\cdot)\) es la función de distribución acumulada de una normal estándar.

3. De manera similar: \[ \textsf{P}(y_{i,j} = 0 \mid \mu, \delta_i, \delta_j) = 1 - \Phi(\mu + \delta_i + \delta_j). \]

Este enfoque linealiza el modelo al sustituir la función de enlace \(\Phi(\cdot)\) por variables auxiliares \(z_{i,j,k}\) con distribución Normal estándar. Esta reformulación simplifica el cálculo de las distribuciones condicionales completas, permitiendo que adopten formas estándar.

Distribuciones Condicionales Completas

1. \(z_{i,j} \mid y_{i,j}, \mu, \delta_i, \delta_j\) es Normal truncada: \[ z_{i,j} \mid y_{i,j}, \mu, \delta_i, \delta_j \sim \begin{cases} \textsf{N}_{(0,\infty)}\,\,\, (\mu + \delta_i + \delta_j, 1), & \text{ si } y_{i,j} = 1, \\ \textsf{N}_{(-\infty,0]}(\mu + \delta_i + \delta_j, 1), & \text{ si } y_{i,j} = 0. \end{cases} \]

2. \(\mu \mid \mathbf{z}, \boldsymbol{\delta}, \sigma^2 \sim \textsf{N}(m, v^2)\), con: \[ v^2 = \left(\frac{1}{\sigma^2} + \sum_{i<j} 1\right)^{-1}, \quad m = v^2 \sum_{i<j} (z_{i,j} - \delta_i - \delta_j). \]

3. \(\delta_i \mid \mathbf{z}, \mu, \tau^2 \sim \textsf{N}(m_i, v_i^2)\), con:
   \[ v_i^2 = \left(\frac{1}{\tau^2} + \sum_{j \neq i} 1\right)^{-1}, \quad m_i = v_i^2 \sum_{j \neq i} (z_{i,j} - \mu - \delta_j). \]

4. \(\sigma^2 \mid \mu \sim \textsf{GI}(a, b)\), con:
   \[ a = a_\sigma + \frac{1}{2}, \quad b = b_\sigma + \frac{\mu^2}{2}. \]

5. \(\tau^2 \mid \boldsymbol{\delta} \sim \textsf{GI}(a, b)\), con:
   \[ a = a_\tau + \frac{n}{2}, \quad b_\tau^* = b_\tau + \frac{\sum_{i} \delta_i^2}{2}. \]

# Distribuciones condicionales completas
# DCC 1: Muestreo de z
sample_z <- function(y, mu, delta, z) {
  for (i in 1:(n - 1)) {
    for (j in (i + 1):n) {
      mean_z <- mu + delta[i] + delta[j]
      if (y[i, j] == 1) {
        z[i, j] <- truncnorm::rtruncnorm(n = 1, a = 0, b = Inf, mean = mean_z, sd = 1)
      } else {
        z[i, j] <- truncnorm::rtruncnorm(n = 1, a = -Inf, b = 0, mean = mean_z, sd = 1)
      }
      z[j, i] <- z[i, j]  # Simetría
    }
  }
  return(z)
}

# DCC 2: Muestreo de mu
sample_mu <- function(z, delta, sigma2) {
  v2_mu <- 1 / (1 / sigma2 + sum(upper.tri(z)))
  m_mu <- v2_mu * sum(z[upper.tri(z)] - delta[row(z)[upper.tri(z)]] - delta[col(z)[upper.tri(z)]])
  return(rnorm(1, mean = m_mu, sd = sqrt(v2_mu)))
}

# DCC 3: Muestreo de delta
sample_delta <- function(z, mu, tau2, delta) {
  for (i in 1:n) {
    neighbors <- setdiff(1:n, i)
    v2_delta <- 1 / (1 / tau2 + length(neighbors))
    m_delta <- v2_delta * sum(z[i, neighbors] - mu - delta[neighbors])
    delta[i] <- rnorm(1, mean = m_delta, sd = sqrt(v2_delta))
  }
  return(delta)
}

# DCC 4: Muestreo de sigma^2
sample_sigma2 <- function(mu) {
  a_sigma_post <- a_sigma + 0.5
  b_sigma_post <- b_sigma + 0.5 * mu^2
  return(1 / rgamma(1, shape = a_sigma_post, rate = b_sigma_post))
}

# DCC 5: Muestreo de tau^2
sample_tau2 <- function(delta) {
  a_tau_post <- a_tau + n / 2
  b_tau_post <- b_tau + 0.5 * sum(delta^2)
  return(1 / rgamma(1, shape = a_tau_post, rate = b_tau_post))
}

# Muestreador de Gibbs
gibbs_sampler <- function(y, n_iter, n_burn, n_thin) {
  # Inicialización
  mu <- 0
  delta <- rnorm(n, 0, 1)
  sigma2 <- 1
  tau2 <- 1
  z <- matrix(0, n, n)  # Variables auxiliares
  # Almacenamiento
  n_samples <- (n_iter - n_burn) / n_thin
  samples <- list(mu = numeric(n_samples), 
                  delta = matrix(0, nrow = n_samples, ncol = n), 
                  sigma2 = numeric(n_samples), 
                  tau2 = numeric(n_samples))
  # Muestreo
  cat("Iniciando muestreador de Gibbs...\n")
  for (t in 1:n_iter) {
    # Llamar a las funciones de muestreo
    z <- sample_z(y, mu, delta, z)
    mu <- sample_mu(z, delta, sigma2)
    delta <- sample_delta(z, mu, tau2, delta)
    sigma2 <- sample_sigma2(mu)
    tau2 <- sample_tau2(delta)
    # Almacenar muestras según n_thin
    if (t > n_burn && (t - n_burn) %% n_thin == 0) {
      idx <- (t - n_burn) / n_thin
      samples$mu[idx] <- mu
      samples$delta[idx, ] <- delta
      samples$sigma2[idx] <- sigma2
      samples$tau2[idx] <- tau2
    }
    # Mostrar progreso
    if (t %% (n_iter / 10) == 0) {
      cat(sprintf("Progreso: %d%% completado\n", (t / n_iter) * 100))
    }
  }
  cat("Muestreador completado.\n")
  return(samples)
}

Hiperparámetros

Se recomienda \(\mathbb{E}[\sigma^2] \leq \mathbb{E}[\tau^2]\) para que la heterogeneidad individual \(\delta_i\) sea mayor que la incertidumbre global \(\mu\).

\(a_\sigma = 2, b_\sigma = 1\).

• Valor esperado: \(\textsf{E}[\sigma^2] = \frac{b_\sigma}{a_\sigma - 1} = 1\).
  
• Coeficiente de variación: \(\textsf{CV}[\sigma^2] = \sqrt{\frac{1}{a_\sigma - 2}} = \infty\).

\(a_\sigma = 1, b_\sigma = 1\).

• Valor esperado* \(\textsf{E}[\sigma^2] = \frac{b_\sigma}{a_\sigma - 1} = \infty\).
  
• Coeficiente de variación: No definido (alta dispersión y cola pesada).

\(a_\tau = 3, b_\tau = 2\).

• Valor esperado: \(\textsf{E}[\tau^2] = \frac{b_\tau}{a_\tau - 1} = 1\).
  
• Coeficiente de variación: \(\textsf{CV}[\tau^2] = \sqrt{\frac{1}{a_\tau - 2}} = \sqrt{\frac{1}{1}} = 1\).

\(a_\tau = 2, b_\tau = 1\).

• Valor esperado: \(\mathbb{E}[\tau^2] = \frac{b_\tau}{a_\tau - 1} = 1\).
  
• Coeficiente de variación: \(\textsf{CV}[\tau^2] = \sqrt{\frac{1}{a_\tau - 2}} = \infty\).

Ejemplo: Simulación

Generación de datos

# Simulación de datos según el modelo
simulate_data <- function(n, mu_true, tau2_true) {
  # Generar parámetros verdaderos
  delta_true <- rnorm(n, mean = 0, sd = sqrt(tau2_true))  # Efectos nodales
  z_true <- matrix(0, n, n)  # Variables latentes
  Y <- matrix(0, n, n)  # Matriz de adyacencia
  # Generar la matriz de adyacencia
  for (i in 1:(n - 1)) {
    for (j in (i + 1):n) {
      z_true[i, j] <- rnorm(1, mean = mu_true + delta_true[i] + delta_true[j], sd = 1)
      z_true[j, i] <- z_true[i, j]  # Simetría
      Y[i, j] <- ifelse(z_true[i, j] > 0, 1, 0)
      Y[j, i] <- Y[i, j]  # Simetría
    }
  }
  # Retornar los datos simulados
  list(Y = Y, z_true = z_true, delta_true = delta_true, mu_true = mu_true, tau2_true = tau2_true)
}

# Parámetros verdaderos
n <- 50
mu_true <- -0.5
tau2_true <- 0.5
# Simulación
set.seed(123)
sim_data <- simulate_data(n, mu_true, tau2_true)
# Matriz de adyacencia simulada
Y <- sim_data$Y

Ajuste del modelo

# Hiperparámetros
a_sigma <- 2 
b_sigma <- 1
a_tau   <- 2 
b_tau   <- 1

# Ajustar el modelo usando Gibbs
n_iter <- 100000 + 10000
n_burn <- 10000
n_thin <- 5
set.seed(123)
samples <- gibbs_sampler(Y, n_iter, n_burn, n_thin)
save(samples, file = "samples_regresion_binaria_sinteticos.RData")

Convergencia

Para cada muestra de \(\mu\) y \(\boldsymbol{\delta}\), la log-verosimilitud del modelo se calcula como:
\[ \log p(\mu, \boldsymbol{\delta} \mid \mathbf{y}) = \sum_{i < j} \left[ y_{i,j} \log(\theta_{i,j}) + (1 - y_{i,j}) \log(1 - \theta_{i,j}) \right], \] donde \(\theta_{i,j} = \Phi(\mu + \delta_i + \delta_j)\).

# Función para calcular la log-verosimilitud
log_likelihood <- function(y, samples) {
  n <- nrow(y)  # Número de nodos
  log_lik_samples <- numeric(length(samples$mu))  # Almacenar la log-verosimilitud para cada muestra
  for (s in seq_along(samples$mu)) {
    mu <- samples$mu[s]
    delta <- samples$delta[s, ]
    log_lik <- 0
    for (i in 1:(n - 1)) {
      for (j in (i + 1):n) {
        eta_ij <- mu + delta[i] + delta[j]  # Predictor lineal
        p_ij <- pnorm(eta_ij)  # Probabilidad del modelo (probit)
        # Sumar la contribución del par (i, j) a la log-verosimilitud
        log_lik <- log_lik + y[i, j] * log(p_ij + 1e-10) + (1 - y[i, j]) * log(1 - p_ij + 1e-10)
      }
    }
    log_lik_samples[s] <- log_lik
  }
  return(log_lik_samples)
}

# Calcular la log-verosimilitud para las muestras del muestreador
log_lik <- log_likelihood(Y, samples)

# Cargar librerías necesarias
suppressMessages(suppressWarnings(library(coda)))

# Función para calcular MCSE
mc_ee <- function(samples) {
  n_eff <- effectiveSize(samples)          # Tamaño efectivo
  sd_est <- sd(as.vector(samples))         # Desviación estándar
  mcse <- abs(sd_est / sqrt(n_eff))        # Error estándar de Monte Carlo (MCSE)
  return(mcse)
}

round(mc_ee(samples$mu), 4)

##   var1 
## 0.0044

round(mc_ee(samples$sigma2), 4)

##   var1 
## 0.0069

round(summary(mc_ee(samples$delta)), 4)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.0068  0.0072  0.0075  0.0075  0.0077  0.0081

Comparación Relativa con la Desviación Estándar:

El MCSE se considera pequeño si es menor o igual al 5% de la desviación estándar posterior (\(\hat{\sigma}\)):
\[ \frac{\text{MCSE}}{\hat{\sigma}} \leq 0.05 \quad (5\%). \]

Comparación con la Media Estimada:

El MCSE se considera pequeño si es menor o igual al 1% de la media estimada (\(\hat{\theta}\)):
\[ \frac{\text{MCSE}}{|\hat{\theta}|} \leq 0.15 \quad (15\%). \]

Evaluación Basada en el Tamaño Efectivo de Muestras (\(n_{\text{eff}}\)):

** El MCSE es inversamente proporcional al tamaño efectivo de muestra (\(n_{\text{eff}}\)):
\[ \text{MCSE} \propto \frac{1}{\sqrt{n_{\text{eff}}}}. \]
Criterio: - \(n_{\text{eff}} > 1000\): MCSE pequeño (alta precisión).
- \(n_{\text{eff}} < 200\): MCSE grande (baja precisión).

Inferencia

# inferencia sobre mu, sigma2 y tau2
# mu_true = -0.5
# tau2_true = 0.5
round(quantile(samples$mu,     probs = c(0.025,0.5,0.975)), 3)

##   2.5%    50%  97.5% 
## -0.742 -0.378 -0.028

round(quantile(samples$tau2,   probs = c(0.025,0.5,0.975)), 3)

##  2.5%   50% 97.5% 
## 0.270 0.413 0.651

# Valores reales
delta_true <- sim_data$delta_true
# Crear un DataFrame con las estimaciones y los valores reales
delta_mean <- colMeans(samples$delta)  # Media posterior de delta
delta_ci95 <- apply(samples$delta, 2, quantile, probs = c(0.025, 0.975))  # IC al 95%
delta_ci99 <- apply(samples$delta, 2, quantile, probs = c(0.005, 0.995))  # IC al 99%
delta_df <- data.frame(
  Node = 1:n,
  Delta_Est = delta_mean,
  Delta_True = delta_true,
  CI95_Lower = delta_ci95[1, ],
  CI95_Upper = delta_ci95[2, ],
  CI99_Lower = delta_ci99[1, ],
  CI99_Upper = delta_ci99[2, ]
)
# Graficar usando ggplot2
ggplot(delta_df, aes(x = Node)) +
  # IC al 99%
  geom_linerange(aes(ymin = CI99_Lower, ymax = CI99_Upper), color = "gray70", linewidth = 1) +
  # IC al 95%
  geom_linerange(aes(ymin = CI95_Lower, ymax = CI95_Upper), color = "lightblue", linewidth = 1.5) +
  # Estimaciones puntuales
  geom_point(aes(y = Delta_Est), color = "blue", size = 3) +
  # Valores reales
  geom_point(aes(y = Delta_True), color = "red", shape = 4, size = 3) +
  # Personalización
  labs(
    title = "Estimaciones de Delta con Intervalos de Credibilidad",
    x = "Vértice",
    y = expression(delta)
  ) +
  theme_minimal(base_size = 14) +
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold")
  )

Ejemplo: Trabajo colaborativo

Los datos de Lazega corresponden a una red de relaciones de trabajo colaborativo entre miembros de una firma de abogados.

Estos datos se recopilaron para estudiar la cooperación entre actores sociales dentro de una organización, analizando el intercambio de diversos tipos de recursos entre ellos.

# librerías
suppressMessages(suppressWarnings(library(igraph)))
suppressMessages(suppressWarnings(library(ggraph)))
suppressMessages(suppressWarnings(library(sand)))
# grafo
g <- graph_from_data_frame(d = elist.lazega, directed = "F")
# clase de objeto
class(g)

## [1] "igraph"

# dirigida?
is_directed(g)

## [1] FALSE

# ponderada?
is_weighted(g)

## [1] FALSE

# orden
(n <- vcount(g))

## [1] 34

# tamaño
(s <- ecount(g))

## [1] 115

# matriz de adyacencia
Y <- as.matrix(as_adjacency_matrix(graph = g, names = F))
# clase de objeto
class(Y)

## [1] "matrix" "array"

# dimensión
dim(Y)

## [1] 34 34

# simétrica?
isSymmetric(Y)

## [1] TRUE

Ajuste del modelo

# Hiperparámetros
a_sigma <- 2 
b_sigma <- 1
a_tau   <- 2 
b_tau   <- 1

# Ajustar el modelo usando Gibbs
n_iter <- 100000 + 10000
n_burn <- 10000
n_thin <- 5
samples <- gibbs_sampler(Y, n_iter, n_burn, n_thin)
save(samples, file = "samples_regresion_binaria_lazega.RData")

# Cargar librerías necesarias
suppressMessages(suppressWarnings(library(coda)))

# Función para calcular MCSE
mc_ee <- function(samples) {
  n_eff <- effectiveSize(samples)          # Tamaño efectivo
  sd_est <- sd(as.vector(samples))         # Desviación estándar
  mcse <- abs(sd_est / sqrt(n_eff))        # Error estándar de Monte Carlo (MCSE)
  return(mcse)
}

round(mc_ee(samples$mu), 4)

##   var1 
## 0.0032

round(mc_ee(samples$sigma2), 4)

##   var1 
## 0.0098

round(summary(mc_ee(samples$delta)), 4)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.0040  0.0042  0.0043  0.0043  0.0044  0.0045

Inferencia

Los parámetros \(\delta_i\) representan la socialidad individual de cada individuo, es decir, su propensión a formar conexiones en la red.

• Valores positivos indican una mayor tendencia a conectarse (alta socialidad).
  
• Valores negativos reflejan una menor tendencia (baja socialidad).
  
• Valores cercanos a cero sugieren un comportamiento promedio.

La varianza \(\tau^2\) mide la heterogeneidad entre los \(\delta_i\):
- Alta varianza \(\tau^2\) implica mayor dispersión en la socialidad (nodos centrales y periféricos).
- Baja varianza \(\tau^2\) sugiere una red socialmente homogénea.

Probabilidades de interacción

La matriz estimada de probabilidades de interacción se calcula iterando sobre las muestras posteriores del modelo y evaluando \(\theta_{ij} = \Phi(\mu + \delta_i + \delta_j)\) para cada par de nodos. Finalmente, se promedian estas matrices a lo largo de todas las muestras para obtener la probabilidad esperada de conexión entre los nodos.

# Función para calcular la matriz de probabilidades theta_ij
compute_theta <- function(samples) {
  n_samples <- length(samples$mu)
  n <- ncol(samples$delta)
  theta_avg <- matrix(0, n, n) # Inicializar matriz promedio
  # Iterar sobre cada muestra
  for (s in 1:n_samples) {
    mu <- samples$mu[s]
    delta <- samples$delta[s, ]
    theta <- matrix(0, n, n)
    # Calcular theta_ij para cada par (i, j)
    for (i in 1:(n - 1)) {
      for (j in (i + 1):n) {
        theta[i, j] <- pnorm(mu + delta[i] + delta[j]) # Probit link
        theta[j, i] <- theta[i, j]  # Simetría
      }
    }
    # Promediar sobre las muestras
    theta_avg <- theta_avg + theta/n_samples
  }
  return(theta_avg)
}

# Calcular la matriz promedio de theta_ij
theta_avg <- compute_theta(samples)

Agrupamiento

Las probabilidades de coagrupamiento (Co-clustering probabilities) se estiman a partir del agrupamiento de los parámetros de sociabilidad \(\delta_i\) en cada iteración del muestreador de Gibbs. El procedimiento utiliza el algoritmo k-means para asignar grupos en cada iteración y calcula las probabilidades de pertenencia conjunta al mismo grupo.

# Librerías
suppressMessages(suppressWarnings(library(cluster)))    # Para la silueta
suppressMessages(suppressWarnings(library(factoextra))) # Para validación de clusters
# Función para determinar automáticamente el número óptimo de clusters
find_optimal_k <- function(delta) {
  max_k <- min(10, length(delta)) # Máximo número de clusters a evaluar
  wss <- numeric(max_k) # Suma de cuadrados dentro de los clusters
  # Calcular WSS para cada k
  for (k in 2:max_k) {
    model <- kmeans(delta, centers = k, nstart = 10)
    wss[k] <- model$tot.withinss # Suma de cuadrados dentro del cluster
  }
  # Método de codo para encontrar el punto óptimo
  optimal_k <- which(diff(diff(wss)) > 0)[1] + 1
  if (is.na(optimal_k)) optimal_k <- 2 # Por defecto al menos 2 clusters
  return(optimal_k)
}
# Función para calcular probabilidades de coagrupamiento
compute_coclustering <- function(samples) {
  n_samples <- length(samples$mu)  # Número de muestras
  n <- ncol(samples$delta)         # Número de nodos
  coclustering <- matrix(0, n, n)  # Inicializar matriz de coagrupamiento
  # Iterar sobre cada muestra
  for (s in 1:n_samples) {
    # Obtener los delta para la muestra actual
    delta <- samples$delta[s, ]
    # Determinar automáticamente el número óptimo de clusters
    k <- find_optimal_k(delta)
    # Aplicar k-means con el k óptimo
    clusters <- kmeans(delta, centers = k, nstart = 10)$cluster
    # Actualizar matriz de co-clustering
    for (i in 1:(n - 1)) {
      for (j in (i + 1):n) {
        if (clusters[i] == clusters[j]) {
          coclustering[i, j] <- coclustering[i, j] + 1
          coclustering[j, i] <- coclustering[i, j] # Simetría
        }
      }
    }
  }
  # Promediar sobre el número total de muestras
  coclustering <- coclustering / n_samples
  return(coclustering)
}

# Calcular probabilidades de coagrupamiento
coclustering_probs <- compute_coclustering(samples)

## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

Bondad de ajuste

La validación del modelo mediante la distribución predictiva posterior evalúa su capacidad para reproducir las características estructurales de la red real.

Se generan redes simuladas en cada iteración del muestreador Gibbs y se calculan estadísticas de prueba como densidad, transitividad, asortatividad, distancia geodésica promedio, grado promedio y desviación estándar del grado.

Estas se comparan con las distribuciones posteriores simuladas, reportando la media posterior y el intervalo de credibilidad al 95%, proporcionando evidencia de ajuste y adecuación del modelo propuesto.

# Calcular las estadísticas observadas en la red real
obs_density <- edge_density(g)  # Densidad
obs_transitivity <- transitivity(g, type = "global")  # Transitividad
obs_assortativity <- assortativity_degree(g, directed = FALSE)  # Asortatividad
obs_avg_path <- mean_distance(g, directed = FALSE)  # Distancia geodésica promedio
obs_avg_degree <- mean(degree(g))  # Grado promedio
obs_sd_degree <- sd(degree(g))  # Desviación estándar del grado
# Guardar en un vector
obs_stats <- c(
  Density = obs_density,
  Transitivity = obs_transitivity,
  Assortativity = obs_assortativity,
  AvgPathLength = obs_avg_path,
  AvgDegree = obs_avg_degree,
  SD_Degree = obs_sd_degree
)

# Inicializar listas para almacenar estadísticas simuladas
sim_density <- c()
sim_transitivity <- c()
sim_assortativity <- c()
sim_avg_path <- c()
sim_avg_degree <- c()
sim_sd_degree <- c()
# Iterar sobre las muestras del muestreador
for (iter in 1:length(samples$mu)) {
  # Extraer parámetros para esta iteración
  mu <- samples$mu[iter]
  delta <- samples$delta[iter, ]
  # Construir la matriz de probabilidad
  P <- pnorm(mu + outer(delta, delta, "+"))
  # Simular la red binaria
  Y_sim <- matrix(rbinom(length(P), size = 1, prob = P), n, n)
  Y_sim[lower.tri(Y_sim)] <- t(Y_sim)[lower.tri(Y_sim)]  # Simetrizar
  diag(Y_sim) <- 0  # Sin bucles
  # Crear grafo a partir de la matriz simulada
  g_sim <- graph_from_adjacency_matrix(Y_sim, mode = "undirected")
  # Calcular estadísticas para la red simulada
  sim_density <- c(sim_density, edge_density(g_sim))
  sim_transitivity <- c(sim_transitivity, transitivity(g_sim, type = "global"))
  sim_assortativity <- c(sim_assortativity, assortativity_degree(g_sim, directed = FALSE))
  sim_avg_path <- c(sim_avg_path, mean_distance(g_sim, directed = FALSE))
  sim_avg_degree <- c(sim_avg_degree, mean(degree(g_sim)))
  sim_sd_degree <- c(sim_sd_degree, sd(degree(g_sim)))
}

## Warning: The dot-dot notation (`..density..`) was deprecated in ggplot2 3.4.0.
## ℹ Please use `after_stat(density)` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

Validación cruzada

La validación cruzada implementada divide los elementos de la matriz de adyacencia \(Y\) (triangular superior) en dos subconjuntos: entrenamiento (80%) y validación (20%). Los valores faltantes en el conjunto de entrenamiento se tratan como datos latentes y se actualizan iterativamente mediante un muestreador de Gibbs, donde en cada iteración se muestrean según la distribución muestral para mantener la coherencia del modelo.

Las observaciones faltantes muestreadas se almacenan para calcular la probabilidad posterior de conexión promediando los valores generados a lo largo del algoritmo. Finalmente, estas probabilidades predictivas se comparan con los datos reales de validación, evaluando el rendimiento del modelo mediante la curva ROC y el AUC para medir su capacidad discriminativa.

# Librerías necesarias
suppressWarnings(
  suppressMessages({
    library(igraph)
    library(ggraph)
    library(tidygraph)
    library(pROC)
    library(ggplot2)
  })
)

# Funciones auxiliares para muestreo Gibbs
gibbs_sampler_missing <- function(Y, n_iter, n_burn, n_thin, valid_indices, indices) {
  n <- nrow(Y)
  # Inicialización de parámetros
  mu <- 0
  delta <- rnorm(n, 0, 1)
  sigma2 <- 1
  tau2 <- 1
  z <- matrix(0, n, n)  # Variables auxiliares
  # Almacenamiento solo para Y_missing
  n_samples <- (n_iter - n_burn) / n_thin
  Y_missing_samples <- matrix(0, nrow = n_samples, ncol = length(valid_indices))
  cat("Iniciando muestreador de Gibbs...
")
  for (t in 1:n_iter) {
    # Paso 1: Muestrear los Y faltantes
    for (idx in valid_indices) {
      i <- indices[idx, 1]
      j <- indices[idx, 2]
      eta_ij <- mu + delta[i] + delta[j]
      prob_ij <- pnorm(eta_ij)
      Y[i, j] <- rbinom(1, 1, prob_ij)
      Y[j, i] <- Y[i, j]  # Simetría
    }
    # Paso 2: Actualizar parámetros condicionales completas
    z <- sample_z(Y, mu, delta, z)
    mu <- sample_mu(z, delta, sigma2)
    delta <- sample_delta(z, mu, tau2, delta)
    sigma2 <- sample_sigma2(mu)
    tau2 <- sample_tau2(delta)
    # Paso 3: Almacenar solo muestras de Y_missing según n_thin
    if (t > n_burn && (t - n_burn) %% n_thin == 0) {
      idx_sample <- (t - n_burn) / n_thin
      Y_missing_samples[idx_sample, ] <- sapply(valid_indices, function(k) {
        i <- indices[k, 1]
        j <- indices[k, 2]
        Y[i, j]
      })
    }
    # Mostrar progreso
    if (t %% (n_iter / 10) == 0) {
      cat(sprintf("Progreso: %d%% completado\n", (t / n_iter) * 100))
    }
  }
  cat("Muestreador completado.\n")
  return(Y_missing_samples)
}

# Hiperparámetros
a_sigma <- 2 
b_sigma <- 1
a_tau   <- 2 
b_tau   <- 1

# Crear folds excluyentes y exhaustivos
set.seed(123)
n <- nrow(Y)
indices <- which(upper.tri(Y, diag = FALSE), arr.ind = TRUE)
n_edges <- nrow(indices)
L <- 5
folds <- split(sample(1:n_edges), rep(1:L, length.out = n_edges))

# Validación cruzada con 5 folds
auc_values <- numeric(L)
roc_list <- list()
for (l in 1:L) {
  cat("Procesando fold", l, "de", L, "\n")
  valid_indices <- folds[[l]]
  Y_train <- Y
  for (idx in valid_indices) {
    i <- indices[idx, 1]
    j <- indices[idx, 2]
    Y_train[i, j] <- NA
    Y_train[j, i] <- NA
  }
  # Muestreador Gibbs
  n_iter <- 100000 + 10000 
  n_burn <- 10000
  n_thin <- 5
  Y_missing_samples <- gibbs_sampler_missing(Y_train, n_iter, n_burn, n_thin, valid_indices, indices)
  # Evaluación de resultados
  Y_valid_values <- Y[indices[valid_indices, ]]
  prob_missing <- colMeans(Y_missing_samples)
  roc_curve <- roc(Y_valid_values, prob_missing)
  auc_values[l] <- auc(roc_curve)
  roc_data <- data.frame(
    FPR = 1 - roc_curve$specificities,
    TPR = roc_curve$sensitivities
  )
  roc_list[[l]] <- roc_data
}
# Guardar resultados
save(auc_values, file = "auc_values.RData")
save(roc_list, file = "roc_list.RData")

# Promedio y desviación estándar del AUC
auc_mean <- mean(auc_values)
auc_sd <- sd(auc_values)
cat("AUC promedio:", round(auc_mean, 3), "\n")

## AUC promedio: 0.714

cat("AUC CV:", round(auc_sd/auc_mean, 3), "\n")

## AUC CV: 0.09

# Graficar todas las curvas ROC en un solo panel
ggplot() +
  geom_line(data = do.call(rbind, lapply(1:L, function(l) cbind(roc_list[[l]], fold = l))),
            aes(x = FPR, y = TPR, color = as.factor(fold)), size = 1, alpha = 0.6) +
  geom_abline(slope = 1, intercept = 0, linetype = "dashed", color = "gray") +
  labs(title = "Curvas ROC para los 5 folds", x = "FPR", y = "TPR", color = "Fold") +
  theme_minimal(base_size = 14)

Modelo con covariables nodales

El modelo para la probabilidad de conexión entre dos nodos se define como:
\[ y_{i,j} \mid \theta_{i,j} \overset{\text{iid}}{\sim} \textsf{Ber}(\theta_{i,j}), \quad \text{para } i < j, \] donde la probabilidad \(\theta_{i,j}\) está dada por:
\[ \theta_{i,j} = \Phi(\mu + \delta_i + \delta_j + \boldsymbol{x}_{i,j}^{\top} \boldsymbol{\beta}), \] donde:

• \(\boldsymbol{x}_{i,j}\): Es un vector de covariables nodales asociado con el par \((i, j)\).
  
• \(\boldsymbol{\beta}\): Representa un vector de coeficientes que cuantifican los efectos de las covariables sobre la probabilidad de conexión.

Las covariables \(\boldsymbol{x}_{i,j}\) deben ser escaladas para garantizar estabilidad numérica y evitar problemas de convergencia durante la estimación. El modelo puede incorporar covariables específicas para cada nodo (\(\boldsymbol{x}_i\) y \(\boldsymbol{x}_j\)) y modelarlas mediante combinaciones lineales, como: \[ \theta_{i,j} = \Phi\left( \mu + \delta_i + \delta_j + \boldsymbol{x}_i^{\top} \boldsymbol{\beta}_1 + \boldsymbol{x}_j^{\top} \boldsymbol{\beta}_2 \right). \] Esto permite diferenciar los efectos de las covariables de los nodos origen y destino en la probabilidad de conexión.