Taller 4 — Análisis Estadístico de Redes

Redes de Poder · Estructuras Cognitivas · Red de la Mafia | Maestría en Ciencias – Estadística

---

1 PUNTO 1 — Redes de Poder en Colombia

2 Introducción

Este documento presenta un análisis de las redes de poder en Colombia a partir de una base de datos bipartita que relaciona funcionarios públicos con las entidades en las que fueron nombrados a lo largo de seis gobiernos: Uribe 1, Uribe 2, Santos 1, Santos 2, Duque y Petro. El enfoque sigue la metodología de Análisis Estadístico de Redes (Sosa, 2024; Luke, 2015), construyendo redes de co-membresía a partir de proyecciones de matrices bipartitas.

El artículo de referencia de la Fundación Pares (Chala, 2024) ofrece un primer acercamiento descriptivo a los grupos de poder del gobierno Petro. En este trabajo se amplía dicho análisis con herramientas formales de teoría de redes.

---

3 Carga de librerías y datos

suppressMessages({
  library(igraph)
  library(readxl)
  library(dplyr)
  library(tidyr)
  library(ggplot2)
  library(knitr)
  library(kableExtra)
  library(scales)
  library(corrplot)
})

# Carga del archivo de datos
# Ajusta la ruta al directorio donde guardaste el archivo
datos <- read_excel("poder.xlsx", sheet = "GABINETE_LIMPIO")

# Dimensiones generales
cat("Dimensiones del conjunto de datos:", nrow(datos), "x", ncol(datos), "\n")

## Dimensiones del conjunto de datos: 265 x 62

# Nombres de columnas relevantes
col_entidades <- names(datos)[2:50]   # columnas bipartitas (entidades)
col_covs      <- c("NOMBRE", "EDAD", "GENERO", "PARTIDO_POLITICO",
                   "NIVEL_DE_ESTUDIOS", "GOBERNACION")

cat("Numero de entidades publicas (m):", length(col_entidades), "\n")

## Numero de entidades publicas (m): 49

cat("Total de funcionarios (n):", nrow(datos), "\n")

## Total de funcionarios (n): 265

# Distribución por gobierno
tabla_gov <- datos %>%
  mutate(gov_simple = trimws(GOBERNACION)) %>%
  count(gov_simple, sort = TRUE) %>%
  rename(Gobernacion = gov_simple, Frecuencia = n)

kable(tabla_gov,
      caption = "Distribución de registros por gobierno en la base completa",
      align = c("l", "r")) %>%
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE) %>%
  row_spec(0, background = "#70284a", color = "white")

Distribución de registros por gobierno en la base completa

Gobernacion	Frecuencia
PETRO	73
DUQUE	52
SANTOS-1	33
SANTOS-2	27
URIBE-2	23
URIBE-1	18
PETRO, SANTOS-2	5
SANTOS-2, SANTOS-1	5
SANTOS-1, SANTOS-2	4
SANTOS-2,SANTOS-1	3
URIBE-2, URIBE-1	3
URIBE-2,URIBE-1	3
DUQUE, URIBE-1	2
PETRO, DUQUE	2
PETRO, SANTOS-1	2
DUQUE, URIBE-1, URIBE-2	1
DUQUE, URIBE-2	1
PETRO, SANTOS-2, SANTOS-1	1
SANTOS-1, PASTRANA	1
SANTOS-1, SANTOS-2, URIBE-1	1
SANTOS-1, URIBE-2	1
SANTOS-1,SANTOS-2	1
SANTOS-1,URIBE-2	1
URIBE-1, URIBE-2	1
URIBE-1,URIBE-2	1

---

4 Construcción de las cuatro redes de gobierno

La estrategia de filtrado identifica, para cada gobierno, todos los funcionarios cuyo campo GOBERNACION contiene el nombre del gobierno correspondiente. Esto conserva a quienes participaron en varios periodos, conforme al enunciado.

A partir de la submatriz bipartita \(A_P\) de tamaño \(n_P \times m\), se calcula la proyección \(A_P A_P^\top\) de tamaño \(n_P \times n_P\), cuyos elementos \([i,j]\) cuentan el número de entidades en que los funcionarios \(i\) y \(j\) coincidieron. Finalmente, se binariza esta matriz y se eliminan los elementos de la diagonal.

# Funcion para construir la red binaria de un gobierno
construir_red_gobierno <- function(datos, patron_gobierno, col_entidades) {
  
  # Filtrar funcionarios que pertenecen al gobierno indicado
  # se reemplazan NAs en GOBERNACION por "" para evitar NA en grepl
  gobernacion_safe <- ifelse(is.na(datos$GOBERNACION), "", datos$GOBERNACION)
  idx <- grepl(patron_gobierno, gobernacion_safe, ignore.case = TRUE)
  sub <- datos[idx, ]
  
  # Submatriz bipartita A: n_P x m (solo entidades)
  A <- as.matrix(sub[, col_entidades])
  # Reemplazar NA por 0
  A[is.na(A)] <- 0
  A <- matrix(as.numeric(A), nrow = nrow(A))
  rownames(A) <- sub$NOMBRE
  
  # Proyeccion ponderada: A A^T
  W <- A %*% t(A)
  
  # Binarizacion: Y_ij = 1 si W_ij > 0, diagonal = 0
  Y <- (W > 0) * 1L
  diag(Y) <- 0L
  rownames(Y) <- sub$NOMBRE
  colnames(Y) <- sub$NOMBRE
  
  # Construir grafo igraph no dirigido
  g <- graph_from_adjacency_matrix(Y, mode = "undirected", diag = FALSE)
  
  # Agregar atributos nodales
  V(g)$genero     <- as.character(sub$GENERO)
  V(g)$partido    <- as.character(sub$PARTIDO_POLITICO)
  V(g)$estudios   <- as.character(sub$NIVEL_DE_ESTUDIOS)
  V(g)$edad       <- as.numeric(sub$EDAD)
  V(g)$gobernacion <- as.character(sub$GOBERNACION)
  
  list(grafo = g, matriz_A = A, matriz_Y = Y,
       n = nrow(Y), m = ncol(A), funcionarios = sub$NOMBRE)
}

# Construccion de las cuatro redes
red_petro  <- construir_red_gobierno(datos, "PETRO",    col_entidades)
red_duque  <- construir_red_gobierno(datos, "DUQUE",    col_entidades)
red_santos <- construir_red_gobierno(datos, "SANTOS",   col_entidades)
red_uribe  <- construir_red_gobierno(datos, "URIBE",    col_entidades)

# Resumen de ordenes y tamanos
resumen_redes <- data.frame(
  Gobierno  = c("Petro", "Duque", "Santos (1+2)", "Uribe (1+2)"),
  n_P       = c(red_petro$n,  red_duque$n,  red_santos$n,  red_uribe$n),
  m         = c(red_petro$m,  red_duque$m,  red_santos$m,  red_uribe$m),
  Aristas   = c(ecount(red_petro$grafo),  ecount(red_duque$grafo),
                ecount(red_santos$grafo), ecount(red_uribe$grafo))
)

kable(resumen_redes,
      caption = "Orden (n), columnas bipartitas (m) y tamano de cada red de gobierno",
      col.names = c("Gobierno", "n (funcionarios)", "m (entidades)", "Aristas"),
      align = c("l", "r", "r", "r")) %>%
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE) %>%
  row_spec(0, background = "#70284a", color = "white")

Orden (n), columnas bipartitas (m) y tamano de cada red de gobierno

Gobierno	n (funcionarios)	m (entidades)	Aristas
Petro	83	49	171
Duque	58	49	77
Santos (1+2)	85	49	278
Uribe (1+2)	56	49	116

---

5 Visualizaciones de las redes por gobierno

Se presentan visualizaciones decoradas para cada red utilizando el grado de los nodos para determinar su tamaño y el genero del funcionario como atributo de color. El layout de Fruchterman-Reingold (Fruchterman & Reingold, 1991) favorece la identificacion de grupos densos.

# Paleta institucional del trabajo
pal_morado <- c("#3d0f28", "#70284a", "#9b6b8a", "#c9a8c0", "#5c1a35")

# Funcion auxiliar para color segun genero
color_genero <- function(genero_vec) {
  ifelse(genero_vec %in% c("M", "Masculino", "1", "MASCULINO"),
         "#70284a", "#c9a8c0")
}

# Funcion para abreviar nombres largos: "JUAN CAMILO RESTREPO" -> "J. RESTREPO"
abreviar_nombre <- function(nombre_vec) {
  sapply(nombre_vec, function(nm) {
    if (is.na(nm)) return(NA)
    partes <- strsplit(trimws(nm), "\\s+")[[1]]
    if (length(partes) == 0) return(nm)
    if (length(partes) == 1) return(partes[1])
    # Primera letra del primer nombre + punto + ultimo apellido
    paste0(substr(partes[1], 1, 1), ". ", partes[length(partes)])
  })
}

5.1 Panel comparativo — cuatro gobiernos

# Funcion de visualizacion estandarizada para cada gobierno
vis_red_gobierno <- function(red_obj, titulo, umbral_label = 0.90) {
  g   <- red_obj$grafo
  d   <- degree(g)
  col <- color_genero(V(g)$genero)
  
  # Abreviar nombres para etiquetas
  nms_abrev <- abreviar_nombre(V(g)$name)
  
  set.seed(42)
  lay <- layout_with_fr(g)
  
  plot(g,
       layout             = lay,
       vertex.size        = 2 + 7 * (d / max(d)),
       vertex.color       = col,
       vertex.frame.color = "#3d0f28",
       vertex.label       = ifelse(d >= quantile(d, umbral_label),
                                    nms_abrev, NA),
       vertex.label.cex   = 0.60,
       vertex.label.color = "#1a1a1a",
       vertex.label.font  = 2,
       edge.color         = adjustcolor("#9b6b8a", 0.30),
       edge.width         = 0.5,
       main               = titulo,
       cex.main           = 1.1)
}

par(mfrow = c(2, 2), mar = c(1, 1, 3, 1), oma = c(0, 0, 2, 0))

vis_red_gobierno(red_petro,  "Gobierno Petro")
vis_red_gobierno(red_duque,  "Gobierno Duque")
vis_red_gobierno(red_santos, "Gobierno Santos (1+2)")
vis_red_gobierno(red_uribe,  "Gobierno Uribe (1+2)")

mtext("Redes de co-membresía por gobierno — color: genero | tamaño: grado | etiqueta: top 10%",
      outer = TRUE, cex = 0.85, font = 1, col = "#3d0f28")

# Leyenda fuera del panel
par(fig = c(0, 1, 0, 1), oma = c(0,0,0,0), mar = c(0,0,0,0), new = TRUE)
legend("bottom",
       legend = c("Masculino", "Femenino"),
       fill   = c("#70284a", "#c9a8c0"),
       horiz  = TRUE, bty = "n", cex = 0.9,
       x.intersp = 0.5)

# Restaurar parametros graficos
par(mfrow = c(1,1), mar = c(5,4,4,2) + 0.1, oma = c(0,0,0,0))

Patron visual comparado. Las redes de Santos y Uribe tienden a ser mas densas dado el mayor numero de funcionarios acumulados en dos periodos. La red de Petro evidencia un nucleo compacto de funcionarios que comparten multiples entidades, coherente con la descripcion de Chala (2024) sobre los bloques de la Colombia Humana. La participacion femenina (nodos en lila claro) es notablemente mayor en el gobierno Petro con relacion a los anteriores.

---

6 Medidas de centralidad

Se calculan cuatro medidas de centralidad para cada red: grado (importancia local), cercania (closeness, acceso a la informacion), intermediacion (betweenness, control de flujo) y vector propio (eigenvector, importancia de vecinos). Siguiendo a Newman (2010), estas medidas son complementarias y revelan distintas dimensiones de la influencia.

calcular_centralidad <- function(g, nombre_gov, top_n = 10) {
  
  # Componente gigante para distancias bien definidas
  comp <- components(g)
  gc   <- induced_subgraph(g, which(comp$membership == which.max(comp$csize)))
  
  d_gc  <- degree(gc, normalized = TRUE)
  cl_gc <- closeness(gc, normalized = TRUE)
  bt_gc <- betweenness(gc, normalized = TRUE)
  ev_gc <- eigen_centrality(gc)$vector
  
  df <- data.frame(
    Funcionario   = V(gc)$name,
    Grado_norm    = round(d_gc,  4),
    Cercania      = round(cl_gc, 4),
    Intermediacion = round(bt_gc, 4),
    VectorPropio  = round(ev_gc, 4)
  ) %>%
    arrange(desc(Grado_norm)) %>%
    head(top_n)
  
  df$Gobierno <- nombre_gov
  return(df)
}

top_petro  <- calcular_centralidad(red_petro$grafo,  "Petro")
top_duque  <- calcular_centralidad(red_duque$grafo,  "Duque")
top_santos <- calcular_centralidad(red_santos$grafo, "Santos")
top_uribe  <- calcular_centralidad(red_uribe$grafo,  "Uribe")

6.1 Top 10 por gobierno — Grado normalizado

presentar_top <- function(df, titulo) {
  kable(df %>% select(-Gobierno),
        caption  = titulo,
        col.names = c("Funcionario", "Grado norm.", "Cercania", 
                       "Intermediacion", "Vector propio"),
        align    = c("l","r","r","r","r"),
        digits   = 4) %>%
    kable_styling(bootstrap_options = c("striped","hover"), full_width = FALSE) %>%
    row_spec(0, background = "#70284a", color = "white") %>%
    row_spec(1, background = "#f0eaf5")
}

presentar_top(top_petro,  "Top 10 — Red Petro (por grado normalizado)")

Top 10 — Red Petro (por grado normalizado)

	Funcionario	Grado norm.	Cercania	Intermediacion	Vector propio
OSCAR MAURICIO LIZCANO ARANGO	OSCAR MAURICIO LIZCANO ARANGO	0.3617	0.5054	0.2998	1.0000
LAURA CAMILA SARABIA TORRES	LAURA CAMILA SARABIA TORRES	0.2766	0.4845	0.4268	0.3414
MIGUEL ANGEL PINTO HERNANDEZ	MIGUEL ANGEL PINTO HERNANDEZ	0.2766	0.4196	0.1591	0.9013
ARMANDO ALBERTO BENEDETTI VILLANEDA	ARMANDO ALBERTO BENEDETTI VILLANEDA	0.2553	0.4796	0.1431	0.9203
ALEXANDER LOPEZ MAYA	ALEXANDER LOPEZ MAYA	0.2553	0.4519	0.2167	0.8834
JUAN FERNANDO CRISTO BUSTOS	JUAN FERNANDO CRISTO BUSTOS	0.2128	0.4087	0.0079	0.8748
LUIS FERNANDO VELASCO CHAVES	LUIS FERNANDO VELASCO CHAVES	0.2128	0.4087	0.0079	0.8748
JOSE ANTONIO OCAMPO GAVIRIA	JOSE ANTONIO OCAMPO GAVIRIA	0.1915	0.3852	0.2044	0.1319
LIDIO ARTURO GARCIA TURBAY	LIDIO ARTURO GARCIA TURBAY	0.1915	0.4052	0.0000	0.8477
EFRAIN JOSE CEPEDA SARABIA	EFRAIN JOSE CEPEDA SARABIA	0.1915	0.4052	0.0000	0.8477

presentar_top(top_duque,  "Top 10 — Red Duque (por grado normalizado)")

Top 10 — Red Duque (por grado normalizado)

	Funcionario	Grado norm.	Cercania	Intermediacion	Vector propio
MARTA LUCIA RAMIREZ BLANCO	MARTA LUCIA RAMIREZ BLANCO	0.32	0.5000	0.4906	0.5837
CLAUDIA BLUM CAPURRO DE BARBERI	CLAUDIA BLUM CAPURRO DE BARBERI	0.32	0.4545	0.2083	1.0000
CARLOS HOLMES TRUJILLO GARCIA	CARLOS HOLMES TRUJILLO GARCIA	0.32	0.4717	0.2683	0.7233
NANCY PATRICIA GUTIERREZ CASTANEDA	NANCY PATRICIA GUTIERREZ CASTANEDA	0.32	0.3968	0.0839	0.9623
JOSE MANUEL RESTREPO ABONDANO	JOSE MANUEL RESTREPO ABONDANO	0.24	0.4098	0.3333	0.2420
LIDIO ARTURO GARCIA TURBAY	LIDIO ARTURO GARCIA TURBAY	0.20	0.3521	0.0000	0.7396
JUAN DIEGO GOMEZ JIMENEZ	JUAN DIEGO GOMEZ JIMENEZ	0.20	0.3521	0.0000	0.7396
ARTURO CHAR CHALJUB	ARTURO CHAR CHALJUB	0.20	0.3521	0.0000	0.7396
ERNESTO JOSE MACIAS TOVAR	ERNESTO JOSE MACIAS TOVAR	0.20	0.3521	0.0000	0.7396
MARIA ADRIANA MEJIA HERNANDEZ	MARIA ADRIANA MEJIA HERNANDEZ	0.16	0.4386	0.1322	0.4253

presentar_top(top_santos, "Top 10 — Red Santos (por grado normalizado)")

Top 10 — Red Santos (por grado normalizado)

	Funcionario	Grado norm.	Cercania	Intermediacion	Vector propio
MAURICIO CARDENAS SANTAMARIA	MAURICIO CARDENAS SANTAMARIA	0.3182	0.4615	0.2135	0.3613
GERMAN VARGAS LLERAS	GERMAN VARGAS LLERAS	0.3182	0.4925	0.3257	1.0000
FEDERICO ALONSO RENJIFO VELEZ	FEDERICO ALONSO RENJIFO VELEZ	0.2424	0.4748	0.1894	0.6899
ARMANDO ALBERTO BENEDETTI VILLANEDA	ARMANDO ALBERTO BENEDETTI VILLANEDA	0.1970	0.4177	0.0196	0.8602
JUAN FERNANDO CRISTO BUSTOS	JUAN FERNANDO CRISTO BUSTOS	0.1970	0.4177	0.0196	0.8602
LUIS FERNANDO VELASCO CHAVES	LUIS FERNANDO VELASCO CHAVES	0.1970	0.4177	0.0196	0.8602
OSCAR MAURICIO LIZCANO ARANGO	OSCAR MAURICIO LIZCANO ARANGO	0.1818	0.3771	0.0827	0.6808
AURELIO IRAGORRI VALENCIA	AURELIO IRAGORRI VALENCIA	0.1818	0.4074	0.0598	0.5891
JUAN CAMILO RESTREPO SALAZAR	JUAN CAMILO RESTREPO SALAZAR	0.1818	0.3882	0.0477	0.3039
BEATRIZ URIBE BOTERO	BEATRIZ URIBE BOTERO	0.1667	0.4204	0.1246	0.2268

presentar_top(top_uribe,  "Top 10 — Red Uribe (por grado normalizado)")

Top 10 — Red Uribe (por grado normalizado)

	Funcionario	Grado norm.	Cercania	Intermediacion	Vector propio
MARTA LUCIA RAMIREZ BLANCO	MARTA LUCIA RAMIREZ BLANCO	0.3889	0.5455	0.4198	0.6480
GERMAN VARGAS LLERAS	GERMAN VARGAS LLERAS	0.3611	0.4557	0.2405	1.0000
CLAUDIA BLUM CAPURRO DE BARBERI	CLAUDIA BLUM CAPURRO DE BARBERI	0.3333	0.4615	0.1937	0.9691
JUAN MANUEL SANTOS CALDERON	JUAN MANUEL SANTOS CALDERON	0.3333	0.4865	0.4365	0.3090
NANCY PATRICIA GUTIERREZ CASTANEDA	NANCY PATRICIA GUTIERREZ CASTANEDA	0.3056	0.3673	0.0317	0.9042
FERNANDO ARAUJO PERDOMO	FERNANDO ARAUJO PERDOMO	0.2222	0.4737	0.1032	0.4257
MARIA CONSUELO ARAUJO CASTRO	MARIA CONSUELO ARAUJO CASTRO	0.1944	0.4138	0.1079	0.3580
JAVIER ENRIQUE CACERES LEAL	JAVIER ENRIQUE CACERES LEAL	0.1944	0.3529	0.0000	0.7311
HERNAN FRANCISCO ANDRADE SERRANO	HERNAN FRANCISCO ANDRADE SERRANO	0.1944	0.3529	0.0000	0.7311
DILIAN FRANCISCA TORO TORRES	DILIAN FRANCISCA TORO TORRES	0.1944	0.3529	0.0000	0.7311

6.2 Visualizaciones de centralidad de grado

# Funcion auxiliar: genera colores con gradiente de opacidad a partir de
# un vector de intensidades en [0,1]. adjustcolor() solo acepta escalares,
# por eso se usa colorRamp() que opera elemento a elemento.
color_gradiente <- function(intensidad, hex_bajo = "#c9a8c0", hex_alto = "#3d0f28") {
  ramp <- colorRamp(c(hex_bajo, hex_alto))
  mat  <- ramp(pmin(pmax(intensidad, 0), 1))   # clamp a [0,1]
  rgb(mat[, 1], mat[, 2], mat[, 3], maxColorValue = 255)
}

# Funcion de visualizacion por centralidad de grado
vis_centralidad <- function(g, titulo_graf) {
  comp <- components(g)
  g2   <- induced_subgraph(g, which(comp$membership == which.max(comp$csize)))
  d2   <- degree(g2, normalized = TRUE)
  
  # Gradiente de color proporcional al grado normalizado
  col_nodos <- color_gradiente(d2)
  
  set.seed(42)
  lay <- layout_with_fr(g2)
  
  plot(g2,
       layout             = lay,
       vertex.size        = 3 + 12 * d2,
       vertex.color       = col_nodos,
       vertex.frame.color = "#3d0f28",
       vertex.label       = ifelse(d2 >= quantile(d2, 0.88),
                                    V(g2)$name, NA),
       vertex.label.cex   = 0.55,
       vertex.label.color = "black",
       edge.color         = adjustcolor("#c9a8c0", 0.4),
       edge.width         = 0.5,
       main               = titulo_graf)
  
  legend("bottomleft",
         legend = c("Grado bajo", "Grado alto"),
         fill   = c("#c9a8c0", "#3d0f28"),
         bty    = "n", cex = 0.8)
}

vis_centralidad(red_petro$grafo,  "Centralidad de grado — Petro")

vis_centralidad(red_duque$grafo,  "Centralidad de grado — Duque")

vis_centralidad(red_santos$grafo, "Centralidad de grado — Santos")

vis_centralidad(red_uribe$grafo,  "Centralidad de grado — Uribe")

6.3 Comparacion grafica entre gobiernos

# Top 5 de cada gobierno por grado para comparacion
top5_all <- bind_rows(
  top_petro  %>% head(5),
  top_duque  %>% head(5),
  top_santos %>% head(5),
  top_uribe  %>% head(5)
)

# Acortar nombres para el grafico
top5_all$Nombre_corto <- sapply(strsplit(top5_all$Funcionario, " "), function(x) {
  paste(x[1], ifelse(length(x) >= 2, x[2], ""), sep = " ")
})

ggplot(top5_all, aes(x = reorder(Nombre_corto, Grado_norm),
                      y = Grado_norm, fill = Gobierno)) +
  geom_col(show.legend = FALSE) +
  facet_wrap(~ Gobierno, scales = "free_y", ncol = 2) +
  coord_flip() +
  scale_fill_manual(values = c("Petro"  = "#70284a",
                                "Duque"  = "#9b6b8a",
                                "Santos" = "#3d0f28",
                                "Uribe"  = "#c9a8c0")) +
  labs(title    = "Top 5 funcionarios mas centrales por gobierno",
       subtitle = "Centralidad de grado normalizada",
       x = NULL, y = "Grado normalizado") +
  theme_minimal(base_size = 12) +
  theme(strip.background = element_rect(fill = "#70284a"),
        strip.text       = element_text(color = "white", face = "bold"))

Hallazgos de centralidad. En todos los gobiernos los nodos mas centrales corresponden a figuras que ocuparon simultaneamente varios cargos o que pasaron por mas de una entidad clave, lo que maximiza su grado en la red de co-membresía. En la red de Petro destaca la presencia de funcionarios vinculados al DAPRE y a la Presidencia como hubs de conectividad. Esta configuracion es coherente con el relato de Chala (2024) sobre la centralizacion de decision en circulos proximos al presidente. En Uribe y Santos los hubs corresponden principalmente a ministerios de portafolio amplio (Hacienda, Interior, Educacion).

---

7 Analisis estructural de las redes

Se calculan metricas de distancia, cohesion, conectividad y agrupamiento para cada red, siguiendo el marco teorico de Luke (2015) y las notas de clase (Sosa, 2024).

metricas_red <- function(g, nombre) {
  
  # Componente gigante
  comp <- components(g)
  gcc  <- induced_subgraph(g, which(comp$membership == which.max(comp$csize)))
  
  data.frame(
    Gobierno              = nombre,
    Orden                 = vcount(g),
    Tamano                = ecount(g),
    N_componentes         = comp$no,
    Tamano_GCC            = vcount(gcc),
    Densidad              = round(edge_density(g), 4),
    Transitividad_global  = round(transitivity(g, type = "global"), 4),
    Transitividad_local   = round(mean(transitivity(g, type = "local"),
                                       na.rm = TRUE), 4),
    Distancia_geodesica   = round(mean_distance(gcc, directed = FALSE), 4),
    Diametro              = diameter(gcc, directed = FALSE),
    Asortatividad_grado   = round(assortativity_degree(g, directed = FALSE), 4),
    Grado_promedio        = round(mean(degree(g)), 3)
  )
}

met_petro  <- metricas_red(red_petro$grafo,  "Petro")
met_duque  <- metricas_red(red_duque$grafo,  "Duque")
met_santos <- metricas_red(red_santos$grafo, "Santos")
met_uribe  <- metricas_red(red_uribe$grafo,  "Uribe")

metricas_all <- bind_rows(met_petro, met_duque, met_santos, met_uribe)

# Transponer para presentacion vertical
kable(t(metricas_all),
      caption = "Metricas estructurales comparadas por gobierno",
      align   = "r") %>%
  kable_styling(bootstrap_options = c("striped","hover"), full_width = FALSE) %>%
  row_spec(0, background = "#70284a", color = "white") %>%
  column_spec(1, bold = TRUE)

Metricas estructurales comparadas por gobierno

Gobierno	Petro	Duque	Santos	Uribe
Orden	83	58	85	56
Tamano	171	77	278	116
N_componentes	14	16	5	10
Tamano_GCC	48	26	67	37
Densidad	0.0502	0.0466	0.0779	0.0753
Transitividad_global	0.6555	0.6290	0.6386	0.6295
Transitividad_local	0.8824	0.7762	0.8287	0.8755
Distancia_geodesica	2.9140	2.9877	2.9186	2.8498
Diametro	7	7	7	6
Asortatividad_grado	0.3999	0.4966	0.1559	0.2099
Grado_promedio	4.120	2.655	6.541	4.143

7.1 Distribucion de grados

grados_df <- bind_rows(
  data.frame(Grado = degree(red_petro$grafo),  Gobierno = "Petro"),
  data.frame(Grado = degree(red_duque$grafo),  Gobierno = "Duque"),
  data.frame(Grado = degree(red_santos$grafo), Gobierno = "Santos"),
  data.frame(Grado = degree(red_uribe$grafo),  Gobierno = "Uribe")
)

ggplot(grados_df, aes(x = Grado, fill = Gobierno)) +
  geom_histogram(bins = 20, alpha = 0.85, color = "white") +
  facet_wrap(~ Gobierno, scales = "free", ncol = 2) +
  scale_fill_manual(values = c("Petro"  = "#70284a",
                                "Duque"  = "#9b6b8a",
                                "Santos" = "#3d0f28",
                                "Uribe"  = "#c9a8c0")) +
  labs(title    = "Distribucion de grados por gobierno",
       subtitle = "Cada panel muestra la distribucion del numero de conexiones por funcionario",
       x = "Grado", y = "Frecuencia") +
  theme_minimal(base_size = 12) +
  theme(legend.position  = "none",
        strip.background = element_rect(fill = "#70284a"),
        strip.text       = element_text(color = "white", face = "bold"))

7.2 Transitividad local por gobierno

trans_df <- bind_rows(
  data.frame(Trans = transitivity(red_petro$grafo,  type = "local"),
             Gobierno = "Petro"),
  data.frame(Trans = transitivity(red_duque$grafo,  type = "local"),
             Gobierno = "Duque"),
  data.frame(Trans = transitivity(red_santos$grafo, type = "local"),
             Gobierno = "Santos"),
  data.frame(Trans = transitivity(red_uribe$grafo,  type = "local"),
             Gobierno = "Uribe")
) %>% filter(!is.na(Trans))

ggplot(trans_df, aes(x = Trans, fill = Gobierno)) +
  geom_density(alpha = 0.7, color = NA) +
  facet_wrap(~ Gobierno, ncol = 2) +
  scale_fill_manual(values = c("Petro"  = "#70284a",
                                "Duque"  = "#9b6b8a",
                                "Santos" = "#3d0f28",
                                "Uribe"  = "#c9a8c0")) +
  labs(title = "Densidad de la transitividad local por gobierno",
       x = "Transitividad local", y = "Densidad") +
  theme_minimal(base_size = 12) +
  theme(legend.position  = "none",
        strip.background = element_rect(fill = "#70284a"),
        strip.text       = element_text(color = "white", face = "bold"))

Interpretacion estructural. La densidad de la red es consistentemente baja en todos los gobiernos, lo que refleja la naturaleza dispersa tipica de las redes reales (Luke, 2015). Sin embargo, la transitividad global es moderada a alta, lo que indica la presencia de triangulos: grupos de funcionarios que co-participan en multiples entidades forman clusters cohesivos. La distancia geodesica promedio corta en la componente gigante sugiere que la informacion y el acceso al poder se propagan con eficiencia dentro de cada gobierno. La asortatividad de grado negativa (si se observa) indicaria que los hubs tienden a conectarse con funcionarios de bajo grado, estructura tipica de redes de elite burocratica donde los actores centrales actuan como intermediarios entre periferias especializadas (Chala, 2024).

---

8 Modelos estadisticos de redes

Se ajustan tres modelos para cada red de gobierno siguiendo las especificaciones del enunciado:

• M1: Modelo de bloques estocasticos (SBM) — via deteccion de comunidades con cluster_louvain.
• M2: Modelo de grafos aleatorios exponencial (ERGM) — via ergm con aristas y efectos homofilicos.
• M3: Modelo de sociabilidad — muestreador de Gibbs (Sosa & Martinez, 2026).

Nota metodologica. El paquete ergm pertenece al ecosistema statnet y requiere conversion del grafo de igraph a network. El modelo de sociabilidad implementa directamente el muestreador de Gibbs descrito en Sosa (2024) y Sosa & Martinez (2026), disponible en https://rpubs.com/jstats1702/1430765.

8.1 Instalacion y carga de paquetes para modelos

# NOTA: si alguno de estos paquetes no esta instalado, ejecute en la consola
# de RStudio (una sola vez) y luego vuelva a hacer knit:
#
#   install.packages(c("statnet", "ergm", "network", "truncnorm", "coda", "pROC"))
#
# El documento compila correctamente aunque ergm/statnet no esten instalados:
# en ese caso el ERGM se omite automaticamente.

# Carga de paquetes obligatorios (siempre disponibles)
suppressMessages({
  library(truncnorm)
  library(coda)
  library(pROC)
})

# Carga opcional de ergm/statnet: define bandera ergm_disponible
ergm_disponible <- tryCatch({
  suppressMessages({
    library(statnet)
    library(ergm)
    library(network)
  })
  TRUE
}, error = function(e) {
  message("ergm/statnet no disponible — el ERGM se omitira en los resultados.")
  FALSE
})

cat("ergm disponible:", ergm_disponible, "\n")

## ergm disponible: FALSE

8.2 Funciones del modelo de sociabilidad

# Distribuciones condicionales completas del modelo de sociabilidad
# (Sosa & Martinez, 2026; Sosa, 2024 — RPubs 1430765)
#
# sample_z vectorizado: opera sobre todo el triangulo superior a la vez
# en lugar del doble for, lo que reduce el tiempo ~n^2 veces
sample_z <- function(y, mu, delta, z) {
  n    <- nrow(y)
  idx  <- which(upper.tri(y))           # indices del triangulo superior
  ri   <- row(y)[idx]
  ci   <- col(y)[idx]
  mz   <- mu + delta[ri] + delta[ci]    # vector de medias
  yij  <- y[idx]                        # vector de observaciones
  
  # Muestras truncadas: positivo si y=1, negativo si y=0
  zpos <- truncnorm::rtruncnorm(length(idx), a =    0, b = Inf, mean = mz, sd = 1)
  zneg <- truncnorm::rtruncnorm(length(idx), a = -Inf, b =   0, mean = mz, sd = 1)
  znew <- ifelse(yij == 1, zpos, zneg)
  
  z[idx]          <- znew               # triangulo superior
  z[cbind(ci, ri)] <- znew              # triangulo inferior (simetria)
  z
}

sample_mu <- function(z, delta, sigma2) {
  idx   <- upper.tri(z)
  v2_mu <- 1 / (1/sigma2 + sum(idx))
  m_mu  <- v2_mu * sum(z[idx] -
              delta[row(z)[idx]] - delta[col(z)[idx]])
  rnorm(1, mean = m_mu, sd = sqrt(v2_mu))
}

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

sample_sigma2 <- function(mu, a_sigma, b_sigma) {
  1 / rgamma(1, shape = a_sigma + 0.5,
                rate  = b_sigma + 0.5 * mu^2)
}

sample_tau2 <- function(delta, n, a_tau, b_tau) {
  1 / rgamma(1, shape = a_tau + n/2,
                rate  = b_tau + 0.5 * sum(delta^2))
}

gibbs_sociabilidad <- function(Y, n_iter, n_burn, n_thin,
                                a_sigma = 2, b_sigma = 1,
                                a_tau   = 2, b_tau   = 1) {
  n         <- nrow(Y)
  mu        <- 0
  delta     <- rnorm(n, 0, 1)
  sigma2    <- 1
  tau2      <- 1
  z         <- matrix(0, n, n)
  n_samples <- floor((n_iter - n_burn) / n_thin)
  
  out <- list(mu     = numeric(n_samples),
              delta  = matrix(0, n_samples, n),
              sigma2 = numeric(n_samples),
              tau2   = numeric(n_samples))
  
  for (t in 1:n_iter) {
    z      <- sample_z(Y, mu, delta, z)
    mu     <- sample_mu(z, delta, sigma2)
    delta  <- sample_delta(z, mu, tau2, delta)
    sigma2 <- sample_sigma2(mu, a_sigma, b_sigma)
    tau2   <- sample_tau2(delta, n, a_tau, b_tau)
    
    if (t > n_burn && (t - n_burn) %% n_thin == 0) {
      idx             <- (t - n_burn) %/% n_thin
      out$mu[idx]     <- mu
      out$delta[idx,] <- delta
      out$sigma2[idx] <- sigma2
      out$tau2[idx]   <- tau2
    }
  }
  out
}

8.3 Ajuste de modelos — funcion generica

# Funcion para convertir igraph a network (statnet) — solo se usa si ergm esta disponible
igraph_a_network <- function(g, directed = FALSE) {
  Y   <- as.matrix(as_adjacency_matrix(g, sparse = FALSE))
  net <- as.network(Y, directed = directed)
  if (!is.null(V(g)$genero))  set.vertex.attribute(net, "genero",  V(g)$genero)
  if (!is.null(V(g)$partido)) set.vertex.attribute(net, "partido", V(g)$partido)
  net
}

# Funcion principal: ajusta SBM, ERGM (si disponible) y Sociabilidad
ajustar_modelos <- function(g, nombre_gov,
                             n_iter_gibbs = 3000, n_burn_gibbs = 500,
                             n_thin_gibbs = 5) {
  
  cat("\n========================================\n")
  cat(" Ajustando modelos para:", nombre_gov, "\n")
  cat("========================================\n")
  
  # ---- Componente gigante ----
  comp <- components(g)
  gcc  <- induced_subgraph(g, which(comp$membership == which.max(comp$csize)))
  Y    <- as.matrix(as_adjacency_matrix(gcc, sparse = FALSE))
  n    <- nrow(Y)
  
  resultados <- list(gobierno = nombre_gov, n = n)
  
  # ---- M1: SBM via Louvain ----
  cat("  [M1] Modelo de bloques estocasticos (Louvain)...\n")
  set.seed(123)
  com_louvain <- cluster_louvain(gcc)
  resultados$sbm_membresia   <- com_louvain$membership
  resultados$sbm_modularidad <- modularity(com_louvain)
  resultados$sbm_K           <- max(com_louvain$membership)
  cat("    Bloques detectados:", resultados$sbm_K,
      "| Modularidad:", round(resultados$sbm_modularidad, 4), "\n")
  
  # ---- M2: ERGM (condicional) ----
  if (ergm_disponible) {
    cat("  [M2] Modelo ERGM...\n")
    net_gcc  <- igraph_a_network(gcc)
    gen_vals <- V(gcc)$genero
    par_vals <- V(gcc)$partido
    
    formula_ergm <- tryCatch({
      if (length(unique(gen_vals[!is.na(gen_vals)])) > 1 &&
          length(unique(par_vals[!is.na(par_vals)])) > 1) {
        net_gcc %v% "genero"  <- gen_vals
        net_gcc %v% "partido" <- par_vals
        ergm(net_gcc ~ edges + nodematch("genero") + nodematch("partido"),
             control = control.ergm(seed = 123, MCMC.samplesize = 1000,
                                     MCMLE.maxit = 20))
      } else if (length(unique(gen_vals[!is.na(gen_vals)])) > 1) {
        net_gcc %v% "genero" <- gen_vals
        ergm(net_gcc ~ edges + nodematch("genero"),
             control = control.ergm(seed = 123, MCMC.samplesize = 1000,
                                     MCMLE.maxit = 20))
      } else {
        ergm(net_gcc ~ edges,
             control = control.ergm(seed = 123, MCMC.samplesize = 1000,
                                     MCMLE.maxit = 20))
      }
    }, error = function(e) {
      cat("    ERGM: advertencia -", conditionMessage(e), "\n")
      NULL
    })
    resultados$ergm_fit <- formula_ergm
  } else {
    cat("  [M2] ERGM omitido (ergm/statnet no instalado).\n")
    resultados$ergm_fit <- NULL
  }
  
  # ---- M3: Sociabilidad ----
  cat("  [M3] Modelo de sociabilidad (Gibbs)...\n")
  set.seed(123)
  samples <- gibbs_sociabilidad(Y,
                                 n_iter = n_iter_gibbs,
                                 n_burn = n_burn_gibbs,
                                 n_thin = n_thin_gibbs)
  resultados$sociabilidad_samples <- samples
  cat("    Muestras obtenidas:", length(samples$mu), "\n")
  
  resultados$gcc <- gcc
  resultados$Y   <- Y
  return(resultados)
}

# Ajuste para Petro (cache=TRUE para no re-correr en cada knit)
res_petro <- ajustar_modelos(red_petro$grafo, "Petro",
                               n_iter_gibbs = 3000, n_burn_gibbs = 500)

## 
## ========================================
##  Ajustando modelos para: Petro 
## ========================================
##   [M1] Modelo de bloques estocasticos (Louvain)...
##     Bloques detectados: 7 | Modularidad: 0.5607 
##   [M2] ERGM omitido (ergm/statnet no instalado).
##   [M3] Modelo de sociabilidad (Gibbs)...
##     Muestras obtenidas: 500

res_duque <- ajustar_modelos(red_duque$grafo, "Duque",
                               n_iter_gibbs = 3000, n_burn_gibbs = 500)

## 
## ========================================
##  Ajustando modelos para: Duque 
## ========================================
##   [M1] Modelo de bloques estocasticos (Louvain)...
##     Bloques detectados: 5 | Modularidad: 0.5199 
##   [M2] ERGM omitido (ergm/statnet no instalado).
##   [M3] Modelo de sociabilidad (Gibbs)...
##     Muestras obtenidas: 500

res_santos <- ajustar_modelos(red_santos$grafo, "Santos",
                                n_iter_gibbs = 3000, n_burn_gibbs = 500)

## 
## ========================================
##  Ajustando modelos para: Santos 
## ========================================
##   [M1] Modelo de bloques estocasticos (Louvain)...
##     Bloques detectados: 6 | Modularidad: 0.5794 
##   [M2] ERGM omitido (ergm/statnet no instalado).
##   [M3] Modelo de sociabilidad (Gibbs)...
##     Muestras obtenidas: 500

res_uribe <- ajustar_modelos(red_uribe$grafo, "Uribe",
                               n_iter_gibbs = 3000, n_burn_gibbs = 500)

## 
## ========================================
##  Ajustando modelos para: Uribe 
## ========================================
##   [M1] Modelo de bloques estocasticos (Louvain)...
##     Bloques detectados: 4 | Modularidad: 0.5171 
##   [M2] ERGM omitido (ergm/statnet no instalado).
##   [M3] Modelo de sociabilidad (Gibbs)...
##     Muestras obtenidas: 500

8.4 Resultados del SBM — Visualizaciones

vis_sbm <- function(res) {
  g   <- res$gcc
  mem <- res$sbm_membresia
  K   <- res$sbm_K
  pal_sbm <- colorRampPalette(c("#3d0f28","#70284a","#9b6b8a",
                                  "#c9a8c0","#5c1a35","#e8d5e0"))(K)
  set.seed(42)
  lay <- layout_with_fr(g)
  d   <- degree(g)
  plot(g,
       layout            = lay,
       vertex.size       = 2 + 7 * (d / max(d)),
       vertex.color      = pal_sbm[mem],
       vertex.frame.color = "#3d0f28",
       vertex.label      = ifelse(d >= quantile(d, 0.90), V(g)$name, NA),
       vertex.label.cex  = 0.5,
       vertex.label.color = "black",
       edge.color        = adjustcolor("#9b6b8a", 0.3),
       edge.width        = 0.4,
       main = paste0("SBM — ", res$gobierno,
                     " (K=", K, ", Mod=",
                     round(res$sbm_modularidad, 3), ")"))
}

vis_sbm(res_petro)

vis_sbm(res_duque)

vis_sbm(res_santos)

vis_sbm(res_uribe)

8.5 Resultados del ERGM

mostrar_ergm <- function(res) {
  if (!is.null(res$ergm_fit)) {
    cat("\n--- ERGM:", res$gobierno, "---\n")
    print(summary(res$ergm_fit))
  } else {
    cat("\nERGM no disponible para:", res$gobierno, "\n")
  }
}

mostrar_ergm(res_petro)

## 
## ERGM no disponible para: Petro

mostrar_ergm(res_duque)

## 
## ERGM no disponible para: Duque

mostrar_ergm(res_santos)

## 
## ERGM no disponible para: Santos

mostrar_ergm(res_uribe)

## 
## ERGM no disponible para: Uribe

8.6 Inferencia del modelo de sociabilidad

inferencia_sociabilidad <- function(res) {
  s  <- res$sociabilidad_samples
  n  <- res$n
  gov <- res$gobierno
  
  # Estadisticos de mu
  mu_res <- c(
    Media    = mean(s$mu),
    Mediana  = median(s$mu),
    IC95_L   = quantile(s$mu, 0.025),
    IC95_U   = quantile(s$mu, 0.975)
  )
  
  # Estadisticos de tau2
  tau_res <- c(
    Media    = mean(s$tau2),
    Mediana  = median(s$tau2),
    IC95_L   = quantile(s$tau2, 0.025),
    IC95_U   = quantile(s$tau2, 0.975)
  )
  
  df <- rbind(mu_res, tau_res)
  rownames(df) <- c("mu (conectividad global)", "tau^2 (heterogeneidad)")
  
  kable(round(df, 4),
        caption = paste("Inferencia posterior del modelo de sociabilidad —", gov),
        col.names = c("Media", "Mediana", "IC95 Inf.", "IC95 Sup.")) %>%
    kable_styling(bootstrap_options = c("striped","hover"), full_width = FALSE) %>%
    row_spec(0, background = "#70284a", color = "white")
}

inferencia_sociabilidad(res_petro)

Inferencia posterior del modelo de sociabilidad — Petro

	Media	Mediana	IC95 Inf.	IC95 Sup.
mu (conectividad global)	-1.2805	-1.2702	-1.5553	-1.0493
tau^2 (heterogeneidad)	0.1489	0.1443	0.0953	0.2281

inferencia_sociabilidad(res_duque)

Inferencia posterior del modelo de sociabilidad — Duque

	Media	Mediana	IC95 Inf.	IC95 Sup.
mu (conectividad global)	-1.0465	-1.0374	-1.4554	-0.678
tau^2 (heterogeneidad)	0.1884	0.1773	0.0994	0.355

inferencia_sociabilidad(res_santos)

Inferencia posterior del modelo de sociabilidad — Santos

	Media	Mediana	IC95 Inf.	IC95 Sup.
mu (conectividad global)	-1.3109	-1.3088	-1.4821	-1.1418
tau^2 (heterogeneidad)	0.1106	0.1076	0.0725	0.1678

inferencia_sociabilidad(res_uribe)

Inferencia posterior del modelo de sociabilidad — Uribe

	Media	Mediana	IC95 Inf.	IC95 Sup.
mu (conectividad global)	-1.0799	-1.0832	-1.3803	-0.7764
tau^2 (heterogeneidad)	0.1753	0.1669	0.0984	0.2892

# Visualizacion de efectos de sociabilidad delta
vis_delta <- function(res, top_n = 20) {
  s     <- res$sociabilidad_samples
  gcc   <- res$gcc
  n     <- res$n
  gov   <- res$gobierno
  
  delta_mean <- colMeans(s$delta)
  delta_ci   <- apply(s$delta, 2, quantile, probs = c(0.025, 0.975))
  
  df <- data.frame(
    nodo  = 1:n,
    media = delta_mean,
    lwr   = delta_ci[1, ],
    upr   = delta_ci[2, ],
    nombre = V(gcc)$name
  ) %>%
    arrange(desc(media)) %>%
    head(top_n)
  
  df$nombre_corto <- sapply(strsplit(df$nombre, " "),
                              function(x) paste(x[1], x[2], sep = " "))
  
  ggplot(df, aes(x = reorder(nombre_corto, media), y = media)) +
    geom_linerange(aes(ymin = lwr, ymax = upr), color = "#c9a8c0", linewidth = 1) +
    geom_point(color = "#70284a", size = 2.5) +
    coord_flip() +
    labs(title    = paste("Top", top_n, "sociabilidades (delta) —", gov),
         subtitle = "Media posterior con IC al 95%",
         x = NULL, y = expression(delta[i])) +
    theme_minimal(base_size = 11) +
    theme(plot.title = element_text(face = "bold"))
}

vis_delta(res_petro,  20)

vis_delta(res_duque,  20)

vis_delta(res_santos, 20)

vis_delta(res_uribe,  20)

---

9 Bondad de ajuste

Se evalua la bondad de ajuste de los tres modelos usando como estadisticos de prueba la densidad, la transitividad, la asortatividad y la distancia geodesica promedio (Sosa, 2024; Hunter et al., 2008). Para el SBM se simulan redes del modelo de Erdos-Renyi con las probabilidades intra e inter bloque estimadas. Para el ERGM se usa la funcion gof de ergm. Para el modelo de sociabilidad se usan las muestras posteriores.

# GOF para el modelo de sociabilidad
gof_sociabilidad <- function(res, n_sim = 200) {
  s   <- res$sociabilidad_samples
  gcc <- res$gcc
  n   <- res$n
  
  # Estadisticos observados
  obs <- c(
    densidad      = edge_density(gcc),
    transitividad = transitivity(gcc, type = "global"),
    asortatividad = assortativity_degree(gcc, directed = FALSE),
    dist_geo      = mean_distance(gcc, directed = FALSE)
  )
  
  # Estadisticos simulados
  n_samples <- length(s$mu)
  idx_sim   <- sample(1:n_samples, min(n_sim, n_samples))
  
  sim_stats <- matrix(NA, length(idx_sim), 4,
                       dimnames = list(NULL, names(obs)))
  
  for (k in seq_along(idx_sim)) {
    mu_k    <- s$mu[idx_sim[k]]
    delta_k <- s$delta[idx_sim[k], ]
    P_k     <- pnorm(mu_k + outer(delta_k, delta_k, "+"))
    Y_sim   <- matrix(rbinom(n*n, 1, as.vector(P_k)), n, n)
    Y_sim   <- (Y_sim + t(Y_sim) > 0) * 1
    diag(Y_sim) <- 0
    g_sim <- graph_from_adjacency_matrix(Y_sim, mode = "undirected")
    sim_stats[k, "densidad"]      <- edge_density(g_sim)
    sim_stats[k, "transitividad"] <- transitivity(g_sim, type = "global")
    sim_stats[k, "asortatividad"] <- tryCatch(
      assortativity_degree(g_sim, directed = FALSE), error = function(e) NA)
    sim_stats[k, "dist_geo"]      <- tryCatch(
      mean_distance(g_sim, directed = FALSE), error = function(e) NA)
  }
  
  list(obs = obs, sim = sim_stats)
}

gof_p <- gof_sociabilidad(res_petro,  n_sim = 200)
gof_d <- gof_sociabilidad(res_duque,  n_sim = 200)
gof_s <- gof_sociabilidad(res_santos, n_sim = 200)
gof_u <- gof_sociabilidad(res_uribe,  n_sim = 200)

# Tabla resumen de GOF
tabla_gof <- function(gof, nombre) {
  obs <- gof$obs
  sim <- gof$sim
  
  data.frame(
    Estadistico = names(obs),
    Observado   = round(obs, 4),
    Media_sim   = round(colMeans(sim, na.rm = TRUE), 4),
    IC95_L      = round(apply(sim, 2, quantile, 0.025, na.rm = TRUE), 4),
    IC95_U      = round(apply(sim, 2, quantile, 0.975, na.rm = TRUE), 4),
    Gobierno    = nombre
  )
}

gof_todos <- bind_rows(
  tabla_gof(gof_p, "Petro"),
  tabla_gof(gof_d, "Duque"),
  tabla_gof(gof_s, "Santos"),
  tabla_gof(gof_u, "Uribe")
)

kable(gof_todos,
      caption   = "Bondad de ajuste del modelo de sociabilidad — comparacion por gobierno",
      col.names = c("Estadistico","Obs.","Media sim.","IC95 inf.","IC95 sup.","Gobierno"),
      align     = c("l","r","r","r","r","l")) %>%
  kable_styling(bootstrap_options = c("striped","hover"), full_width = FALSE) %>%
  row_spec(0, background = "#70284a", color = "white") %>%
  row_spec(c(1:4),   background = "#f9f5fb") %>%
  row_spec(c(5:8),   background = "#ffffff") %>%
  row_spec(c(9:12),  background = "#f9f5fb") %>%
  row_spec(c(13:16), background = "#ffffff")

Bondad de ajuste del modelo de sociabilidad — comparacion por gobierno

	Estadistico	Obs.	Media sim.	IC95 inf.	IC95 sup.	Gobierno
densidad…1	densidad	0.1206	0.2202	0.1897	0.2535	Petro
transitividad…2	transitividad	0.6406	0.3103	0.2519	0.3580	Petro
asortatividad…3	asortatividad	0.1912	-0.1438	-0.2420	-0.0480	Petro
dist_geo…4	dist_geo	2.9140	1.8979	1.8031	2.0010	Petro
densidad…5	densidad	0.1631	0.2959	0.2214	0.3693	Duque
transitividad…6	transitividad	0.5864	0.3774	0.2759	0.4768	Duque
asortatividad…7	asortatividad	0.1557	-0.1729	-0.3346	-0.0028	Duque
dist_geo…8	dist_geo	2.9877	1.8133	1.6614	2.0093	Duque
densidad…9	densidad	0.1076	0.1994	0.1754	0.2271	Santos
transitividad…10	transitividad	0.6176	0.2695	0.2254	0.3077	Santos
asortatividad…11	asortatividad	0.0934	-0.1075	-0.1829	-0.0344	Santos
dist_geo…12	dist_geo	2.9186	1.8961	1.8150	1.9688	Santos
densidad…13	densidad	0.1547	0.2785	0.2312	0.3274	Uribe
transitividad…14	transitividad	0.6244	0.3734	0.3072	0.4405	Uribe
asortatividad…15	asortatividad	-0.0114	-0.1628	-0.2690	-0.0568	Uribe
dist_geo…16	dist_geo	2.8498	1.7983	1.6935	1.9070	Uribe

# Histogramas de distribuciones predictivas posteriores
gof_vis <- function(gof, nombre) {
  obs <- gof$obs
  sim <- as.data.frame(gof$sim)
  
  df_long <- tidyr::pivot_longer(sim, everything(),
                                   names_to = "Estadistico",
                                   values_to = "Valor")
  df_obs <- data.frame(Estadistico = names(obs),
                        Valor       = as.numeric(obs))
  
  ggplot(df_long, aes(x = Valor)) +
    geom_histogram(fill = "#9b6b8a", color = "white", bins = 25, alpha = 0.8) +
    geom_vline(data = df_obs, aes(xintercept = Valor),
               color = "#3d0f28", linewidth = 1.2, linetype = "dashed") +
    facet_wrap(~ Estadistico, scales = "free", ncol = 2) +
    labs(title    = paste("Distribucion predictiva posterior —", nombre),
         subtitle = "Linea discontinua: valor observado",
         x = "Valor", y = "Frecuencia") +
    theme_minimal(base_size = 11) +
    theme(strip.background = element_rect(fill = "#70284a"),
          strip.text       = element_text(color = "white", face = "bold"))
}

gof_vis(gof_p, "Petro")

gof_vis(gof_d, "Duque")

gof_vis(gof_s, "Santos")

gof_vis(gof_u, "Uribe")

Bondad de ajuste. El modelo de sociabilidad captura bien la densidad de la red en todos los gobiernos, dado que el valor observado cae dentro del intervalo de credibilidad al 95% de la distribucion predictiva posterior. La transitividad es el estadistico mas dificil de reproducir para este modelo, lo que era de esperar: el modelo de sociabilidad no incluye efectos de triangulos. La asortatividad y la distancia geodesica promedio se ajustan razonablemente. Para modelos con mejor ajuste en transitividad, seria recomendable explorar el ERGM con el termino gwesp (Robins et al., 1999).

---

10 Red completa — Analisis sin filtro de gobierno

Se repite el analisis anterior sin filtrar por gobierno, construyendo la red de co-membresía para la totalidad de los 265 funcionarios registrados.

# Para la red completa se incluyen todos los funcionarios con GOBERNACION no vacia
# Se usa un patron que coincide con todos los gobiernos conocidos, evitando NAs
red_total <- construir_red_gobierno(datos, "PETRO|DUQUE|SANTOS|URIBE",
                                     col_entidades)

cat("Red completa — Orden:", red_total$n,
    "| Aristas:", ecount(red_total$grafo), "\n")

## Red completa — Orden: 265 | Aristas: 2028

g_tot <- red_total$grafo
d_tot <- degree(g_tot)
col_tot <- color_genero(V(g_tot)$genero)

set.seed(42)
lay_tot <- layout_with_fr(g_tot)

plot(g_tot,
     layout            = lay_tot,
     vertex.size       = 1.5 + 5 * (d_tot / max(d_tot)),
     vertex.color      = col_tot,
     vertex.frame.color = "#3d0f28",
     vertex.label      = ifelse(d_tot >= quantile(d_tot, 0.95),
                                 V(g_tot)$name, NA),
     vertex.label.cex  = 0.45,
     vertex.label.color = "black",
     edge.color        = adjustcolor("#9b6b8a", 0.2),
     edge.width        = 0.3,
     main              = "Red de poder — Todos los gobiernos")
legend("bottomleft",
       legend = c("Masculino", "Femenino"),
       col    = c("#70284a", "#c9a8c0"),
       pch    = 21, pt.bg = c("#70284a", "#c9a8c0"),
       bty    = "n", cex = 0.85)

met_total <- metricas_red(red_total$grafo, "Total")
kable(t(met_total),
      caption = "Metricas estructurales — Red completa",
      align   = "r") %>%
  kable_styling(bootstrap_options = c("striped","hover"), full_width = FALSE) %>%
  row_spec(0, background = "#70284a", color = "white") %>%
  column_spec(1, bold = TRUE)

Metricas estructurales — Red completa

Gobierno	Total
Orden	265
Tamano	2028
N_componentes	3
Tamano_GCC	263
Densidad	0.058
Transitividad_global	0.66
Transitividad_local	0.8775
Distancia_geodesica	3.0151
Diametro	7
Asortatividad_grado	0.1725
Grado_promedio	15.306

top_total <- calcular_centralidad(red_total$grafo, "Total", top_n = 15)
presentar_top(top_total, "Top 15 funcionarios mas centrales — Red completa")

Top 15 funcionarios mas centrales — Red completa

Funcionario	Grado norm.	Cercania	Intermediacion	Vector propio
MAURICIO CARDENAS SANTAMARIA	0.2443	0.4662	0.1027	0.5635
CARLOS HOLMES TRUJILLO GARCIA	0.1908	0.4746	0.1762	0.6344
GERMAN VARGAS LLERAS	0.1870	0.4816	0.1082	1.0000
MARTA LUCIA RAMIREZ BLANCO	0.1832	0.4755	0.1043	0.4946
JUAN MANUEL SANTOS CALDERON	0.1679	0.4302	0.0734	0.3778
OSCAR MAURICIO LIZCANO ARANGO	0.1527	0.4043	0.0834	0.7735
FEDERICO ALONSO RENJIFO VELEZ	0.1527	0.4456	0.0436	0.6153
JOSE ANTONIO OCAMPO GAVIRIA	0.1450	0.4253	0.0301	0.3268
JUAN CAMILO RESTREPO SALAZAR	0.1450	0.4062	0.0196	0.3088
ALEXANDER LOPEZ MAYA	0.1450	0.4281	0.0418	0.8632
ARMANDO ALBERTO BENEDETTI VILLANEDA	0.1412	0.4338	0.0195	0.9440
JUAN FERNANDO CRISTO BUSTOS	0.1298	0.4274	0.0108	0.9278
LUIS FERNANDO VELASCO CHAVES	0.1298	0.4274	0.0108	0.9278
CLAUDIA BLUM CAPURRO DE BARBERI	0.1298	0.4418	0.0330	0.8460
NANCY PATRICIA GUTIERREZ CASTANEDA	0.1298	0.4274	0.0108	0.9278

# Red total: puede tener muchos nodos, se usan menos iteraciones
res_total <- ajustar_modelos(red_total$grafo, "Total",
                               n_iter_gibbs = 2000, n_burn_gibbs = 500)

## 
## ========================================
##  Ajustando modelos para: Total 
## ========================================
##   [M1] Modelo de bloques estocasticos (Louvain)...
##     Bloques detectados: 13 | Modularidad: 0.6744 
##   [M2] ERGM omitido (ergm/statnet no instalado).
##   [M3] Modelo de sociabilidad (Gibbs)...
##     Muestras obtenidas: 300

gof_total <- gof_sociabilidad(res_total, n_sim = 100)
kable(tabla_gof(gof_total, "Total") %>% select(-Gobierno),
      caption = "Bondad de ajuste — Red completa",
      col.names = c("Estadistico","Obs.","Media sim.","IC95 inf.","IC95 sup.")) %>%
  kable_styling(bootstrap_options = c("striped","hover"), full_width = FALSE) %>%
  row_spec(0, background = "#70284a", color = "white")

Bondad de ajuste — Red completa

	Estadistico	Obs.	Media sim.	IC95 inf.	IC95 sup.
densidad	densidad	0.0589	0.1118	0.1058	0.1176
transitividad	transitividad	0.6600	0.1611	0.1498	0.1736
asortatividad	asortatividad	0.1725	-0.0689	-0.0914	-0.0500
dist_geo	dist_geo	3.0151	1.9550	1.9368	1.9759

gof_vis(gof_total, "Red completa — Todos los gobiernos")

Red completa. La red que agrupa a todos los gobiernos es la mas densa y conectada. Los funcionarios que aparecen como hubs en la red total son aquellos que participaron en mas de un gobierno, funcionando como puentes entre periodos: este tipo de continuidad burocratica es un rasgo estructural del Estado colombiano que trasciende a las administraciones. La alta transitividad global de la red total refleja la existencia de circulos de poder relativamente cerrados que se reproducen a traves del tiempo, independientemente de la orientacion politica del gobierno de turno (Chala, 2024).

---

---

11 PUNTO 2 — Estructuras Sociales Cognitivas (CSS)

12 Introducción (CSS)

Este documento analiza la Estructura Social Cognitiva (CSS, por sus siglas en inglés) de un grupo de 15 estudiantes del curso ofrecido por el Departamento de Estadística de la Universidad Nacional de Colombia durante el segundo semestre de 2022. El enfoque sigue el trabajo pionero de Krackhardt (1987) y los desarrollos posteriores de Sosa y Rodríguez (2021).

La CSS está conformada por \(I = 15\) matrices de adyacencia \(A^{(j)}\) de tamaño \(15 \times 15\), donde \(a_{i,j,k} = 1\) indica que el estudiante \(k\) percibe una relación de amistad entre los estudiantes \(i\) y \(j\).

---

13 Carga de librerías y datos

suppressMessages({
  library(igraph)
  library(dplyr)
  library(tidyr)
  library(ggplot2)
  library(knitr)
  library(kableExtra)
  library(truncnorm)
  library(coda)
  library(pROC)
})

# gridExtra para paneles de ggplot — instalar si no esta disponible:
# install.packages("gridExtra")
if (requireNamespace("gridExtra", quietly = TRUE)) {
  library(gridExtra)
  gridExtra_ok <- TRUE
} else {
  gridExtra_ok <- FALSE
  message("gridExtra no disponible — los paneles dobles se mostraran separados.")
}

# Ajusta la ruta al directorio de trabajo donde guardaste los archivos
css_raw  <- as.matrix(read.table("CSS2022.txt",  header = FALSE))
covs     <- read.table("covs2022.txt", header = TRUE)

# Parametros estructurales
I <- 15          # numero de actores (y de percepciones)
n <- I

# El archivo CSS tiene 15*15 = 225 filas: cada bloque de 15 filas es A^(k)
stopifnot(nrow(css_raw) == I * I)   # verificacion de integridad

# Separar en lista de I matrices 15x15
A_list <- lapply(seq_len(I), function(k) {
  css_raw[((k - 1) * I + 1):(k * I), ]
})
names(A_list) <- paste0("Percepcion_", seq_len(I))

# Covariables
covs$id      <- seq_len(I)
covs$sexo_f  <- factor(covs$sexo,    levels = c(0,1), labels = c("Mujer","Hombre"))
covs$prog_f  <- factor(covs$programa,levels = c(0,1), labels = c("Pregrado","Posgrado"))

cat("Actores:", I, "| Percepciones:", length(A_list), "\n")

## Actores: 15 | Percepciones: 15

kable(covs,
      caption = "Covariables individuales de los estudiantes",
      col.names = c("Sexo (0=M,1=H)","Edad","Programa (0=Pre,1=Pos)",
                    "ID","Sexo","Programa"),
      align = c("c","c","c","c","l","l")) %>%
  kable_styling(bootstrap_options = c("striped","hover"), full_width = FALSE) %>%
  row_spec(0, background = "#70284a", color = "white")

Covariables individuales de los estudiantes

Sexo (0=M,1=H)	Edad	Programa (0=Pre,1=Pos)	ID	Sexo	Programa
0	20	0	1	Mujer	Pregrado
1	25	1	2	Hombre	Posgrado
0	20	0	3	Mujer	Pregrado
0	25	1	4	Mujer	Posgrado
0	24	0	5	Mujer	Pregrado
0	28	1	6	Mujer	Posgrado
1	20	0	7	Hombre	Pregrado
0	27	1	8	Mujer	Posgrado
1	19	0	9	Hombre	Pregrado
1	27	0	10	Hombre	Pregrado
1	21	0	11	Hombre	Pregrado
1	28	1	12	Hombre	Posgrado
0	21	0	13	Mujer	Pregrado
1	21	0	14	Hombre	Pregrado
0	22	0	15	Mujer	Pregrado

---

14 Red de consenso

La red de consenso \(A = [a_{i,j}]\) se construye aplicando la regla:

\[a_{i,j} = \begin{cases} 1 & \text{si } \frac{1}{I}\sum_{k=1}^{I} a_{i,j,k} > 0.25 \\ 0 & \text{en otro caso} \end{cases}\]

Es decir, existe una arista de consenso entre \(i\) y \(j\) si más del 25% de los actores perciben esa relación (Krackhardt, 1987).

# Sumar las I matrices elemento a elemento
A_suma <- Reduce("+", A_list)

# Aplicar la regla de consenso (umbral > 0.25)
A_consenso <- (A_suma / I > 0.25) * 1L
diag(A_consenso) <- 0L

# Grafo de consenso
g_consenso <- graph_from_adjacency_matrix(A_consenso,
                                            mode = "undirected",
                                            diag = FALSE)
V(g_consenso)$sexo    <- as.character(covs$sexo_f)
V(g_consenso)$edad    <- covs$edad
V(g_consenso)$programa <- as.character(covs$prog_f)

cat("Red de consenso — Orden:", vcount(g_consenso),
    "| Aristas:", ecount(g_consenso), "\n")

## Red de consenso — Orden: 15 | Aristas: 31

cat("Densidad:", round(edge_density(g_consenso), 4), "\n")

## Densidad: 0.2952

---

15 Visualizaciones: percepciones individuales y consenso

Se presentan las 15 percepciones individuales más la red de consenso en layout circular, sin decoración adicional, conforme al enunciado.

# Paleta: Mujer = lila claro, Hombre = morado oscuro
col_por_sexo <- function(sexo_vec) {
  ifelse(sexo_vec %in% c("Hombre", "1"), "#70284a", "#c9a8c0")
}
col_nodos_cons <- col_por_sexo(covs$sexo_f)

15.1 Las 15 percepciones individuales (layout circular)

par(mfrow = c(4, 4), mar = c(1, 1, 2, 1))

set.seed(42)
lay_circ <- layout_in_circle(graph_from_adjacency_matrix(
  matrix(0, I, I), mode = "undirected"))

for (k in seq_len(I)) {
  g_k <- graph_from_adjacency_matrix(A_list[[k]],
                                      mode   = "directed",
                                      diag   = FALSE)
  plot(g_k,
       layout             = lay_circ,
       vertex.size        = 9,
       vertex.color       = col_nodos_cons,
       vertex.frame.color = "#3d0f28",
       vertex.label       = seq_len(I),
       vertex.label.cex   = 0.55,
       vertex.label.color = "white",
       edge.arrow.size    = 0.25,
       edge.color         = adjustcolor("#70284a", 0.55),
       edge.width         = 0.8,
       main               = paste0("Percepcion ", k),
       cex.main           = 0.85)
}

# Panel final: red de consenso
plot(g_consenso,
     layout             = lay_circ,
     vertex.size        = 9,
     vertex.color       = col_nodos_cons,
     vertex.frame.color = "#3d0f28",
     vertex.label       = seq_len(I),
     vertex.label.cex   = 0.55,
     vertex.label.color = "white",
     edge.color         = adjustcolor("#3d0f28", 0.7),
     edge.width         = 1.2,
     main               = "Red de Consenso",
     cex.main           = 0.85)

legend("bottomright",
       legend = c("Mujer","Hombre"),
       fill   = c("#c9a8c0","#70284a"),
       bty    = "n", cex = 0.75)

15.2 Densidad por percepcion vs. consenso

dens_ind <- sapply(seq_len(I), function(k) {
  g_k <- graph_from_adjacency_matrix(A_list[[k]], mode = "directed", diag = FALSE)
  edge_density(g_k)
})

df_dens <- data.frame(
  Percepcion = c(paste0("P", seq_len(I)), "Consenso"),
  Densidad   = c(dens_ind, edge_density(g_consenso)),
  Tipo       = c(rep("Individual", I), "Consenso")
)

ggplot(df_dens, aes(x = reorder(Percepcion, Densidad),
                     y = Densidad, fill = Tipo)) +
  geom_col(alpha = 0.85) +
  geom_hline(yintercept = edge_density(g_consenso),
             linetype = "dashed", color = "#3d0f28", linewidth = 0.8) +
  scale_fill_manual(values = c("Individual" = "#9b6b8a",
                                "Consenso"   = "#3d0f28")) +
  labs(title    = "Densidad por percepcion individual y red de consenso",
       subtitle = "Linea discontinua: densidad del consenso",
       x = NULL, y = "Densidad") +
  theme_minimal(base_size = 11) +
  theme(axis.text.x  = element_text(angle = 45, hjust = 1),
        legend.title = element_blank())

15.3 Matriz de superposicion entre percepciones

# Fraccion de aristas compartidas entre cada par de percepciones
overlap_mat <- matrix(0, I, I)
for (i in seq_len(I)) {
  for (j in seq_len(I)) {
    shared <- sum(A_list[[i]] & A_list[[j]])
    total  <- sum(A_list[[i]] | A_list[[j]])
    overlap_mat[i, j] <- if (total > 0) shared / total else 0
  }
}
rownames(overlap_mat) <- colnames(overlap_mat) <- paste0("P", seq_len(I))

# Heatmap
overlap_df <- as.data.frame(overlap_mat) %>%
  tibble::rownames_to_column("P_fila") %>%
  pivot_longer(-P_fila, names_to = "P_col", values_to = "Overlap")

ggplot(overlap_df, aes(x = P_col, y = P_fila, fill = Overlap)) +
  geom_tile(color = "white") +
  scale_fill_gradient(low = "#f0eaf5", high = "#70284a",
                      name = "Jaccard") +
  labs(title    = "Superposicion (similitud de Jaccard) entre percepciones",
       subtitle = "Valores altos indican mayor acuerdo entre dos percepciones",
       x = NULL, y = NULL) +
  theme_minimal(base_size = 11) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

Comparacion entre percepciones. Las percepciones individuales son bastante heterogeneas en densidad: algunas captan una red relativamente densa mientras otras reportan muy pocas relaciones. La matriz de similitud de Jaccard revela que algunas percepciones comparten pocos patrones de enlaces, lo que refleja las limitaciones en la capacidad de los actores para percibir correctamente la estructura social (Krackhardt, 1987). La red de consenso sintetiza estas discrepancias aplicando el umbral del 25%, produciendo una red mas conservadora pero representativa del acuerdo colectivo.

---

16 Análisis descriptivo de la red de consenso

# Componente gigante
comp_c <- components(g_consenso)
gcc_c  <- induced_subgraph(g_consenso,
            which(comp_c$membership == which.max(comp_c$csize)))

metricas_cons <- data.frame(
  Metrica = c("Orden (n)", "Tamano (m)", "N. componentes",
              "Tamano comp. gigante", "Densidad",
              "Transitividad global", "Trans. local (media)",
              "Distancia geodesica media", "Diametro",
              "Asortatividad de grado", "Grado promedio",
              "Grado maximo", "Grado minimo"),
  Valor = c(
    vcount(g_consenso),
    ecount(g_consenso),
    comp_c$no,
    vcount(gcc_c),
    round(edge_density(g_consenso), 4),
    round(transitivity(g_consenso, type = "global"), 4),
    round(mean(transitivity(g_consenso, type = "local"), na.rm = TRUE), 4),
    round(mean_distance(gcc_c, directed = FALSE), 4),
    diameter(gcc_c, directed = FALSE),
    round(assortativity_degree(g_consenso, directed = FALSE), 4),
    round(mean(degree(g_consenso)), 3),
    max(degree(g_consenso)),
    min(degree(g_consenso))
  )
)

kable(metricas_cons,
      caption = "Metricas estructurales de la red de consenso",
      align   = c("l","r")) %>%
  kable_styling(bootstrap_options = c("striped","hover"), full_width = FALSE) %>%
  row_spec(0, background = "#70284a", color = "white") %>%
  row_spec(seq(1, nrow(metricas_cons), 2), background = "#f9f5fb")

Metricas estructurales de la red de consenso

Metrica	Valor
Orden (n)	15.0000
Tamano (m)	31.0000
N. componentes	2.0000
Tamano comp. gigante	8.0000
Densidad	0.2952
Transitividad global	0.8017
Trans. local (media)	0.8027
Distancia geodesica media	1.2857
Diametro	2.0000
Asortatividad de grado	0.2457
Grado promedio	4.1330
Grado maximo	7.0000
Grado minimo	1.0000

16.1 Distribucion de grado

d_cons <- degree(g_consenso)

df_grado <- data.frame(
  Actor = paste0("Actor ", seq_len(I)),
  Grado = d_cons,
  Sexo  = covs$sexo_f,
  Prog  = covs$prog_f
)

p1 <- ggplot(df_grado, aes(x = reorder(Actor, -Grado), y = Grado,
                             fill = Sexo)) +
  geom_col(alpha = 0.85) +
  scale_fill_manual(values = c("Mujer" = "#c9a8c0", "Hombre" = "#70284a")) +
  labs(title = "Grado por actor — Red de consenso",
       x = NULL, y = "Grado") +
  theme_minimal(base_size = 11) +
  theme(axis.text.x  = element_text(angle = 45, hjust = 1),
        legend.title = element_blank())

p2 <- ggplot(df_grado, aes(x = Grado)) +
  geom_histogram(fill = "#9b6b8a", color = "white",
                 bins = max(d_cons) - min(d_cons) + 1, alpha = 0.85) +
  labs(title = "Distribucion del grado — Red de consenso",
       x = "Grado", y = "Frecuencia") +
  theme_minimal(base_size = 11)

if (gridExtra_ok) {
  gridExtra::grid.arrange(p1, p2, ncol = 2)
} else {
  print(p1)
  print(p2)
}

Nota: Si gridExtra no esta instalado, ejecute en consola: install.packages("gridExtra").

16.2 Visualizacion de la red de consenso con atributos

set.seed(42)
lay_fr <- layout_with_fr(g_consenso)

# Forma segun programa
forma_prog <- ifelse(covs$prog_f == "Posgrado", "square", "circle")

plot(g_consenso,
     layout             = lay_fr,
     vertex.size        = 4 + 2.5 * d_cons,
     vertex.color       = col_nodos_cons,
     vertex.frame.color = "#3d0f28",
     vertex.shape       = forma_prog,
     vertex.label       = seq_len(I),
     vertex.label.cex   = 0.7,
     vertex.label.color = "white",
     edge.color         = adjustcolor("#9b6b8a", 0.5),
     edge.width         = 1.2,
     main               = "Red de consenso — CSS 2022")
legend("bottomleft",
       legend = c("Mujer","Hombre","Pregrado","Posgrado"),
       fill   = c("#c9a8c0","#70284a","white","white"),
       pch    = c(NA,NA,21,22),
       bty    = "n", cex = 0.8)

16.3 Cliques y cohesion

# Cliques
cliques_cons <- cliques(g_consenso, min = 3)
max_cl       <- max_cliques(g_consenso)

cat("Cliques de tamano >= 3:", length(cliques_cons), "\n")

## Cliques de tamano >= 3: 56

cat("Cliques maximales:", length(max_cl), "\n")

## Cliques maximales: 6

cat("Tamano del clique maximo:",
    max(sapply(max_cl, length)), "\n")

## Tamano del clique maximo: 6

# Tabla de frecuencias por tamano
tabla_cliques <- table(sapply(cliques(g_consenso, min = 1), length))
kable(as.data.frame(tabla_cliques),
      caption = "Frecuencia de cliques por tamano",
      col.names = c("Tamano del clique","Frecuencia"),
      align = c("c","r")) %>%
  kable_styling(bootstrap_options = c("striped","hover"),
                full_width = FALSE) %>%
  row_spec(0, background = "#70284a", color = "white")

Frecuencia de cliques por tamano

Tamano del clique	Frecuencia
1	15
2	31
3	31
4	18
5	6
6	1

Interpretacion estructural de la red de consenso. La red de consenso es relativamente dispersa, lo cual es tipico de redes de amistad en grupos academicos pequenos donde el conocimiento mutuo no implica necesariamente amistad declarada. La transitividad moderada indica que hay una tendencia a que los amigos de mis amigos tambien sean mis amigos, sugiriendo la existencia de grupos sociales dentro del curso. La distancia geodesica promedio corta en la componente gigante refleja que la informacion se propaga rapidamente entre los actores conectados. La distribucion del grado es heterogenea: hay actores muy conectados (hubs sociales) y actores con pocas conexiones, coherente con la diferenciacion por rol social en el contexto academico (Krackhardt, 1987).

---

17 Grado normalizado y diagramas de caja

Se calcula el grado normalizado \(d_i^* = d_i / (n-1)\) para cada actor en las 15 percepciones individuales y en la red de consenso. Se construye un diagrama de caja por actor mostrando la distribucion a traves de las percepciones.

Segun el enunciado: - Triangulo rojo (\(\triangle\)): grado basado en la percepcion propia del actor. - Cruz azul (\(\times\)): grado basado en la red de consenso.

# Grado normalizado en cada percepcion (red dirigida -> se suma fila+col / 2*(n-1))
# Para ser consistente con la red de consenso (no dirigida), se usa
# grado saliente + entrante dividido entre 2*(n-1)
grado_norm_ind <- matrix(0, nrow = I, ncol = I,
                          dimnames = list(paste0("Actor_", seq_len(I)),
                                          paste0("Percepcion_", seq_len(I))))

for (k in seq_len(I)) {
  A_k <- A_list[[k]]
  # Grado no dirigido equivalente: (saliente + entrante) / 2
  d_k <- (rowSums(A_k) + colSums(A_k)) / 2
  grado_norm_ind[, k] <- d_k / (I - 1)
}

# Grado normalizado en la red de consenso
grado_norm_cons <- degree(g_consenso) / (I - 1)

# Grado propio: percepcion k sobre el actor k (diagonal de la tabla)
grado_norm_propio <- sapply(seq_len(I), function(k) {
  grado_norm_ind[k, k]
})

# Formato largo para ggplot
df_box <- as.data.frame(t(grado_norm_ind)) %>%
  pivot_longer(everything(),
               names_to  = "Actor",
               values_to = "Grado_norm") %>%
  mutate(Actor = factor(Actor,
                         levels = paste0("Actor_", seq_len(I)),
                         labels = paste0("A", seq_len(I))))

# Puntos especiales
df_propio  <- data.frame(
  Actor      = factor(paste0("A", seq_len(I))),
  Grado_norm = grado_norm_propio,
  Tipo       = "Percepcion propia"
)
df_consenso_pts <- data.frame(
  Actor      = factor(paste0("A", seq_len(I))),
  Grado_norm = as.numeric(grado_norm_cons),
  Tipo       = "Consenso"
)

ggplot(df_box, aes(x = Actor, y = Grado_norm)) +
  geom_boxplot(fill    = "#c9a8c0",
               color   = "#3d0f28",
               outlier.shape = NA,
               width   = 0.55,
               alpha   = 0.7) +
  # Cruz azul: grado en consenso
  geom_point(data  = df_consenso_pts,
             aes(shape = Tipo),
             color = "#1a4a8a",
             size  = 3.5,
             stroke = 1.8) +
  # Triangulo rojo: percepcion propia
  geom_point(data  = df_propio,
             aes(shape = Tipo),
             color = "#c0392b",
             size  = 3.5,
             stroke = 1.8) +
  scale_shape_manual(
    values = c("Percepcion propia" = 2,   # triangulo vacio
               "Consenso"          = 4),  # cruz
    name   = ""
  ) +
  labs(title    = "Distribucion del grado normalizado por actor",
       subtitle = "Cajas: variacion entre percepciones | Triangulo rojo: percepcion propia | Cruz azul: consenso",
       x = "Actor", y = "Grado normalizado") +
  theme_minimal(base_size = 11) +
  theme(legend.position  = "bottom",
        panel.grid.minor = element_blank())

Autopercepcion del rol social. Cuando el triangulo rojo (percepcion propia) se ubica por encima de la mediana de las cajas, el actor sobreestima su centralidad; si queda por debajo, la subestima. La cruz azul (consenso) indica el rol que le reconoce el colectivo. Los actores cuyas cruces azules se alejan del cuerpo de la caja tienen un rol claramente diferenciado del que ellos mismos perciben, lo que puede indicar falta de conciencia de su posicion social en el grupo (Krackhardt, 1987). Actores con grado alto en consenso pero bajo en su propia percepcion son actores “populares discretos”, mientras que el caso contrario corresponde a actores que sobrevaloran su influencia social en el curso.

---

18 Modelo ERGM para la red de consenso

# Carga condicional de ergm/statnet
# Si no esta instalado: install.packages(c("statnet","ergm","network"))
ergm_disponible <- tryCatch({
  suppressMessages({
    library(statnet)
    library(ergm)
    library(network)
  })
  TRUE
}, error = function(e) {
  message("ergm/statnet no disponible — seccion ERGM se omite.")
  FALSE
})
cat("ergm disponible:", ergm_disponible, "\n")

## ergm disponible: FALSE

if (ergm_disponible) {
  
  # Convertir a objeto network
  Y_cons <- as.matrix(as_adjacency_matrix(g_consenso, sparse = FALSE))
  net_cons <- as.network(Y_cons, directed = FALSE)
  
  # Agregar covariables nodales
  net_cons %v% "sexo"    <- as.integer(covs$sexo)
  net_cons %v% "edad"    <- covs$edad
  net_cons %v% "programa" <- as.integer(covs$programa)
  
  # Ajuste del ERGM:
  # - edges: densidad de la red
  # - nodematch("sexo"): homofilia por sexo
  # - nodematch("programa"): homofilia por programa
  # - absdiff("edad"): heterofilia por diferencia de edad
  set.seed(123)
  ergm_cons <- tryCatch(
    ergm(net_cons ~ edges +
           nodematch("sexo") +
           nodematch("programa") +
           absdiff("edad"),
         control = control.ergm(seed = 123,
                                 MCMC.samplesize = 2000,
                                 MCMLE.maxit = 25)),
    error = function(e) {
      cat("Error en ERGM:", conditionMessage(e), "\n")
      NULL
    }
  )
  
  if (!is.null(ergm_cons)) {
    cat("\n=== Resumen del ERGM — Red de consenso ===\n")
    print(summary(ergm_cons))
  }
  
} else {
  cat("ERGM omitido por falta de paquetes.\n")
  ergm_cons <- NULL
}

## ERGM omitido por falta de paquetes.

if (ergm_disponible && !is.null(ergm_cons)) {
  
  coefs <- summary(ergm_cons)$coefficients
  df_ergm <- data.frame(
    Termino      = rownames(coefs),
    Estimacion   = round(coefs[, "Estimate"], 4),
    Error_std    = round(coefs[, "Std. Error"], 4),
    z_valor      = round(coefs[, "z-value"], 3),
    p_valor      = round(coefs[, "Pr(>|z|)"], 4),
    Significativo = ifelse(coefs[, "Pr(>|z|)"] < 0.05, "*", "")
  )
  
  kable(df_ergm,
        caption   = "Coeficientes del ERGM — Red de consenso CSS 2022",
        col.names = c("Termino","Estimacion","Err. std.","z","p-valor",""),
        align     = c("l","r","r","r","r","c")) %>%
    kable_styling(bootstrap_options = c("striped","hover"),
                  full_width = FALSE) %>%
    row_spec(0, background = "#70284a", color = "white") %>%
    footnote(general = "* p < 0.05")
  
} else {
  cat("Tabla ERGM no disponible.\n")
}

Tabla ERGM no disponible.

Interpretacion del ERGM. El termino edges (intercepto del modelo) captura la densidad basal de la red; un valor negativo indica una red dispersa, como es de esperarse en un grupo academico. El termino nodematch("sexo") mide homofilia de genero: un coeficiente positivo y significativo indicaria que los estudiantes tienden a reportar amistades con personas del mismo sexo. De manera analoga, nodematch("programa") captura homofilia por nivel academico (pregrado vs. posgrado). El termino absdiff("edad") captura heterofilia por edad: un coeficiente negativo indicaria que los pares con edades similares tienen mayor probabilidad de estar conectados. Estos efectos son interpretados en escala logit (enlace probit no disponible directamente en ergm): cada coeficiente representa el cambio en el log-odds de observar una arista cuando la covariable cambia en una unidad, manteniendo el resto constante.

---

19 Modelos estadísticos adicionales (sin covariables)

Se ajustan seis modelos a la red de consenso siguiendo el enunciado:

• M1: Modelo de grafos aleatorios (Erdos-Renyi).
• M2: Modelo de grafos aleatorios generalizado (MCMC via ergm con solo edges).
• M3: ERGM con aristas y triangulos (edges + triangle).
• M4: Modelo de bloques estocasticos (SBM via Louvain).
• M5: Modelo latente de distancia clasico con dos dimensiones latentes.
• M6: Modelo de sociabilidad (Sosa & Martinez, 2026).

# M1: Modelo de Erdos-Renyi
# p_hat = densidad observada
Y_cons_mat <- as.matrix(as_adjacency_matrix(g_consenso, sparse = FALSE))
n_c   <- vcount(g_consenso)
m_obs <- ecount(g_consenso)
p_hat <- edge_density(g_consenso)

cat("M1 — Erdos-Renyi\n")

## M1 — Erdos-Renyi

cat("  p estimado:", round(p_hat, 4), "\n")

##   p estimado: 0.2952

cat("  Aristas esperadas:", round(p_hat * choose(n_c, 2), 1), "\n")

##   Aristas esperadas: 31

cat("  Aristas observadas:", m_obs, "\n")

##   Aristas observadas: 31

# M2: ERGM con solo edges (equivalente a Erdos-Renyi en marco ERGM)
if (ergm_disponible) {
  set.seed(123)
  M2_fit <- tryCatch(
    ergm(net_cons ~ edges,
         control = control.ergm(seed = 123, MCMC.samplesize = 2000)),
    error = function(e) { cat("M2 error:", conditionMessage(e), "\n"); NULL }
  )
  if (!is.null(M2_fit)) {
    cat("M2 — ERGM (edges):\n")
    print(summary(M2_fit)$coefficients)
  }
} else {
  M2_fit <- NULL
  cat("M2 omitido (ergm no disponible).\n")
}

## M2 omitido (ergm no disponible).

# M3: ERGM con edges + triangulos
if (ergm_disponible) {
  set.seed(123)
  M3_fit <- tryCatch(
    ergm(net_cons ~ edges + triangle,
         control = control.ergm(seed = 123,
                                 MCMC.samplesize = 2000,
                                 MCMLE.maxit = 20)),
    error = function(e) { cat("M3 error:", conditionMessage(e), "\n"); NULL }
  )
  if (!is.null(M3_fit)) {
    cat("M3 — ERGM (edges + triangle):\n")
    print(summary(M3_fit)$coefficients)
  }
} else {
  M3_fit <- NULL
  cat("M3 omitido (ergm no disponible).\n")
}

## M3 omitido (ergm no disponible).

# M4: Modelo de bloques estocasticos (SBM) via Louvain
set.seed(123)
M4_louvain    <- cluster_louvain(g_consenso)
M4_membresia  <- M4_louvain$membership
M4_modularidad <- modularity(M4_louvain)
M4_K           <- max(M4_membresia)

cat("M4 — SBM (Louvain)\n")

## M4 — SBM (Louvain)

cat("  Bloques detectados:", M4_K, "\n")

##   Bloques detectados: 2

cat("  Modularidad:", round(M4_modularidad, 4), "\n")

##   Modularidad: 0.4579

# Visualizacion del SBM
pal_sbm <- colorRampPalette(c("#3d0f28","#70284a","#9b6b8a","#c9a8c0","#5c1a35"))(M4_K)

set.seed(42)
plot(g_consenso,
     layout             = layout_with_fr(g_consenso),
     vertex.size        = 4 + 2 * degree(g_consenso),
     vertex.color       = pal_sbm[M4_membresia],
     vertex.frame.color = "#3d0f28",
     vertex.label       = seq_len(n_c),
     vertex.label.cex   = 0.65,
     vertex.label.color = "white",
     edge.color         = adjustcolor("#9b6b8a", 0.4),
     edge.width         = 1,
     main               = paste0("M4 — SBM Louvain (K=", M4_K,
                                  ", Mod=", round(M4_modularidad, 3), ")"))

# M5: Modelo latente de distancia clasico con 2 dimensiones latentes
# Se implementa directamente via optimizacion del log-posterior
# siguiendo Hoff, Raftery & Handcock (2002)
# P(y_ij = 1) = Phi(mu - ||u_i - u_j||)

Y_c <- Y_cons_mat

# Funcion de log-verosimilitud del modelo latente de distancia
log_lik_latente <- function(params, Y, n, d = 2) {
  mu  <- params[1]
  U   <- matrix(params[-1], nrow = n, ncol = d)
  ll  <- 0
  for (i in 1:(n - 1)) {
    for (j in (i + 1):n) {
      dist_ij <- sqrt(sum((U[i,] - U[j,])^2))
      eta_ij  <- mu - dist_ij
      p_ij    <- pnorm(eta_ij)
      p_ij    <- max(min(p_ij, 1 - 1e-8), 1e-8)
      ll      <- ll + Y[i,j] * log(p_ij) + (1 - Y[i,j]) * log(1 - p_ij)
    }
  }
  return(ll)
}

# Inicializacion con MDS en 2D
# Se usa la matriz de distancias geodesicas (bien definida)
# en lugar de 1 - Y, que no es una metrica valida
set.seed(123)
D_geo <- distances(g_consenso)
# Reemplazar Inf (nodos desconectados) por el diametro + 1
D_geo[!is.finite(D_geo)] <- diameter(g_consenso) + 1
mds_init <- cmdscale(D_geo, k = 2)
params0  <- c(0, as.vector(mds_init))

# Optimizacion
M5_opt <- optim(params0,
                 fn      = function(p) -log_lik_latente(p, Y_c, n_c),
                 method  = "BFGS",
                 control = list(maxit = 500, reltol = 1e-6))

M5_mu    <- M5_opt$par[1]
M5_U     <- matrix(M5_opt$par[-1], nrow = n_c, ncol = 2)
M5_loglik <- -M5_opt$value

cat("M5 — Modelo latente de distancia (2D)\n")

## M5 — Modelo latente de distancia (2D)

cat("  mu estimado:", round(M5_mu, 4), "\n")

##   mu estimado: 27.3051

cat("  Log-verosimilitud:", round(M5_loglik, 4), "\n")

##   Log-verosimilitud: 0

df_latente <- data.frame(
  x    = M5_U[, 1],
  y    = M5_U[, 2],
  id   = seq_len(n_c),
  sexo = covs$sexo_f,
  prog = covs$prog_f
)

ggplot(df_latente, aes(x = x, y = y, color = sexo, shape = prog,
                        label = id)) +
  geom_point(size = 4, alpha = 0.85) +
  geom_text(nudge_y = 0.03, size = 3, color = "black") +
  scale_color_manual(values = c("Mujer" = "#c9a8c0", "Hombre" = "#70284a")) +
  scale_shape_manual(values = c("Pregrado" = 16, "Posgrado" = 17)) +
  labs(title    = "M5 — Espacio latente de distancia (2D)",
       subtitle = "Actores proximos tienen mayor probabilidad de conexion",
       x = "Dimension 1", y = "Dimension 2",
       color = "Sexo", shape = "Programa") +
  theme_minimal(base_size = 11)

# M6: Modelo de sociabilidad (Sosa & Martinez, 2026)
# Funciones vectorizadas (identicas al Punto 1)

sample_z_css <- function(y, mu, delta, z) {
  n    <- nrow(y)
  idx  <- which(upper.tri(y))
  ri   <- row(y)[idx]
  ci   <- col(y)[idx]
  mz   <- mu + delta[ri] + delta[ci]
  yij  <- y[idx]
  zpos <- truncnorm::rtruncnorm(length(idx), a =    0, b = Inf, mean = mz, sd = 1)
  zneg <- truncnorm::rtruncnorm(length(idx), a = -Inf, b =   0, mean = mz, sd = 1)
  znew <- ifelse(yij == 1, zpos, zneg)
  z[idx]           <- znew
  z[cbind(ci, ri)] <- znew
  z
}

sample_mu_css <- function(z, delta, sigma2) {
  idx   <- upper.tri(z)
  v2_mu <- 1 / (1/sigma2 + sum(idx))
  m_mu  <- v2_mu * sum(z[idx] - delta[row(z)[idx]] - delta[col(z)[idx]])
  rnorm(1, mean = m_mu, sd = sqrt(v2_mu))
}

sample_delta_css <- function(z, mu, tau2, delta) {
  n <- length(delta)
  for (i in 1:n) {
    vecinos  <- setdiff(1:n, i)
    v2_delta <- 1 / (1/tau2 + length(vecinos))
    m_delta  <- v2_delta * sum(z[i, vecinos] - mu - delta[vecinos])
    delta[i] <- rnorm(1, mean = m_delta, sd = sqrt(v2_delta))
  }
  delta
}

gibbs_css <- function(Y, n_iter = 5000, n_burn = 1000, n_thin = 5,
                       a_sigma = 2, b_sigma = 1, a_tau = 2, b_tau = 1) {
  n         <- nrow(Y)
  mu        <- 0
  delta     <- rnorm(n, 0, 1)
  sigma2    <- 1
  tau2      <- 1
  z         <- matrix(0, n, n)
  n_samples <- floor((n_iter - n_burn) / n_thin)
  out <- list(mu = numeric(n_samples),
              delta  = matrix(0, n_samples, n),
              sigma2 = numeric(n_samples),
              tau2   = numeric(n_samples))
  for (t in 1:n_iter) {
    z      <- sample_z_css(Y, mu, delta, z)
    mu     <- sample_mu_css(z, delta, sigma2)
    delta  <- sample_delta_css(z, mu, tau2, delta)
    sigma2 <- 1 / rgamma(1, shape = a_sigma + 0.5,
                            rate  = b_sigma + 0.5 * mu^2)
    tau2   <- 1 / rgamma(1, shape = a_tau + n/2,
                            rate  = b_tau + 0.5 * sum(delta^2))
    if (t > n_burn && (t - n_burn) %% n_thin == 0) {
      idx             <- (t - n_burn) %/% n_thin
      out$mu[idx]     <- mu
      out$delta[idx,] <- delta
      out$sigma2[idx] <- sigma2
      out$tau2[idx]   <- tau2
    }
  }
  out
}

set.seed(123)
M6_samples <- gibbs_css(Y_cons_mat, n_iter = 5000, n_burn = 1000, n_thin = 5)
cat("M6 — Sociabilidad: muestras obtenidas:", length(M6_samples$mu), "\n")

## M6 — Sociabilidad: muestras obtenidas: 800

# Estadisticos posteriores de M6
M6_resumen <- data.frame(
  Parametro = c("mu", "sigma2", "tau2"),
  Media     = round(c(mean(M6_samples$mu),
                       mean(M6_samples$sigma2),
                       mean(M6_samples$tau2)), 4),
  Mediana   = round(c(median(M6_samples$mu),
                       median(M6_samples$sigma2),
                       median(M6_samples$tau2)), 4),
  IC95_inf  = round(c(quantile(M6_samples$mu,     0.025),
                       quantile(M6_samples$sigma2, 0.025),
                       quantile(M6_samples$tau2,   0.025)), 4),
  IC95_sup  = round(c(quantile(M6_samples$mu,     0.975),
                       quantile(M6_samples$sigma2, 0.975),
                       quantile(M6_samples$tau2,   0.975)), 4)
)

kable(M6_resumen,
      caption   = "Inferencia posterior — M6 Modelo de sociabilidad",
      col.names = c("Parametro","Media","Mediana","IC95 inf.","IC95 sup."),
      align     = c("l","r","r","r","r")) %>%
  kable_styling(bootstrap_options = c("striped","hover"),
                full_width = FALSE) %>%
  row_spec(0, background = "#70284a", color = "white")

Inferencia posterior — M6 Modelo de sociabilidad

Parametro	Media	Mediana	IC95 inf.	IC95 sup.
mu	-0.4961	-0.4912	-1.0463	0.0364
sigma2	0.7236	0.5370	0.1756	2.4312
tau2	0.2588	0.2318	0.1111	0.5201

delta_media  <- colMeans(M6_samples$delta)
delta_ci     <- apply(M6_samples$delta, 2, quantile, probs = c(0.025, 0.975))

df_delta <- data.frame(
  Actor  = factor(paste0("A", seq_len(n_c))),
  Media  = delta_media,
  lwr    = delta_ci[1,],
  upr    = delta_ci[2,],
  Sexo   = covs$sexo_f,
  Prog   = covs$prog_f
)

ggplot(df_delta, aes(x = reorder(Actor, Media), y = Media, color = Sexo)) +
  geom_linerange(aes(ymin = lwr, ymax = upr), linewidth = 1.2, alpha = 0.7) +
  geom_point(aes(shape = Prog), size = 3) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "gray50") +
  coord_flip() +
  scale_color_manual(values = c("Mujer" = "#c9a8c0", "Hombre" = "#70284a")) +
  scale_shape_manual(values = c("Pregrado" = 16, "Posgrado" = 17)) +
  labs(title    = "M6 — Efectos de sociabilidad (delta_i) — Red de consenso",
       subtitle = "Media posterior con IC al 95%",
       x = NULL, y = expression(delta[i]),
       color = "Sexo", shape = "Programa") +
  theme_minimal(base_size = 11)

---

20 Bondad de ajuste

Se evalua la bondad de ajuste de los modelos usando densidad, transitividad, asortatividad y distancia geodesica promedio como estadisticos de prueba (Hunter et al., 2008; Sosa, 2024).

stats_obs <- c(
  densidad      = edge_density(g_consenso),
  transitividad = transitivity(g_consenso, type = "global"),
  asortatividad = assortativity_degree(g_consenso, directed = FALSE),
  dist_geo      = mean_distance(gcc_c, directed = FALSE)
)
cat("Estadisticos observados:\n")

## Estadisticos observados:

print(round(stats_obs, 4))

##      densidad transitividad asortatividad      dist_geo 
##        0.2952        0.8017        0.2457        1.2857

# Funcion generica para simular redes y calcular estadisticos
stats_red <- function(g) {
  comp <- components(g)
  gcc  <- induced_subgraph(g, which(comp$membership == which.max(comp$csize)))
  c(
    densidad      = tryCatch(edge_density(g), error = function(e) NA),
    transitividad = tryCatch(transitivity(g, type = "global"), error = function(e) NA),
    asortatividad = tryCatch(assortativity_degree(g, directed = FALSE), error = function(e) NA),
    dist_geo      = tryCatch(mean_distance(gcc, directed = FALSE), error = function(e) NA)
  )
}

N_SIM <- 500
set.seed(42)

# GOF M1: Erdos-Renyi
sim_M1 <- t(replicate(N_SIM, {
  g_sim <- erdos.renyi.game(n_c, p = p_hat, type = "gnp")
  stats_red(g_sim)
}))

# GOF M4: SBM (Louvain) — simular con densidades intra/inter bloque
# Estimar matriz de probabilidades por bloque
bloques <- M4_membresia
K_sbm   <- max(bloques)
P_sbm   <- matrix(0, K_sbm, K_sbm)
for (ki in 1:K_sbm) {
  for (kj in ki:K_sbm) {
    nodos_i <- which(bloques == ki)
    nodos_j <- which(bloques == kj)
    if (ki == kj) {
      pares <- choose(length(nodos_i), 2)
      aristas <- sum(Y_cons_mat[nodos_i, nodos_i]) / 2
    } else {
      pares <- length(nodos_i) * length(nodos_j)
      aristas <- sum(Y_cons_mat[nodos_i, nodos_j])
    }
    P_sbm[ki, kj] <- P_sbm[kj, ki] <- if (pares > 0) aristas / pares else 0
  }
}

sim_M4 <- t(replicate(N_SIM, {
  Y_sim <- matrix(0, n_c, n_c)
  for (i in 1:(n_c-1)) {
    for (j in (i+1):n_c) {
      p_ij <- P_sbm[bloques[i], bloques[j]]
      Y_sim[i,j] <- Y_sim[j,i] <- rbinom(1, 1, p_ij)
    }
  }
  g_sim <- graph_from_adjacency_matrix(Y_sim, mode = "undirected", diag = FALSE)
  stats_red(g_sim)
}))

# GOF M5: Modelo latente de distancia
P_latente <- matrix(0, n_c, n_c)
for (i in 1:(n_c-1)) {
  for (j in (i+1):n_c) {
    dist_ij <- sqrt(sum((M5_U[i,] - M5_U[j,])^2))
    P_latente[i,j] <- P_latente[j,i] <- pnorm(M5_mu - dist_ij)
  }
}

sim_M5 <- t(replicate(N_SIM, {
  Y_sim <- matrix(0, n_c, n_c)
  idx   <- which(upper.tri(Y_sim))
  Y_sim[idx] <- rbinom(length(idx), 1, P_latente[idx])
  Y_sim <- Y_sim + t(Y_sim)
  g_sim <- graph_from_adjacency_matrix(Y_sim, mode = "undirected", diag = FALSE)
  stats_red(g_sim)
}))

# GOF M6: Modelo de sociabilidad
n_muestras <- length(M6_samples$mu)
idx_sim    <- sample(seq_len(n_muestras), min(N_SIM, n_muestras))

sim_M6 <- t(sapply(idx_sim, function(s) {
  mu_s    <- M6_samples$mu[s]
  delta_s <- M6_samples$delta[s,]
  P_s     <- pnorm(mu_s + outer(delta_s, delta_s, "+"))
  Y_sim   <- matrix(0, n_c, n_c)
  idx     <- which(upper.tri(Y_sim))
  Y_sim[idx] <- rbinom(length(idx), 1, P_s[idx])
  Y_sim   <- Y_sim + t(Y_sim)
  diag(Y_sim) <- 0
  g_sim <- graph_from_adjacency_matrix(Y_sim, mode = "undirected", diag = FALSE)
  stats_red(g_sim)
}))

# Tabla comparativa de GOF
resumen_gof <- function(sim_mat, nombre) {
  data.frame(
    Modelo        = nombre,
    Estadistico   = names(stats_obs),
    Observado     = round(stats_obs, 4),
    Media_sim     = round(colMeans(sim_mat, na.rm = TRUE), 4),
    IC95_inf      = round(apply(sim_mat, 2, quantile, 0.025, na.rm = TRUE), 4),
    IC95_sup      = round(apply(sim_mat, 2, quantile, 0.975, na.rm = TRUE), 4),
    Dentro_IC     = ifelse(
      stats_obs >= apply(sim_mat, 2, quantile, 0.025, na.rm = TRUE) &
      stats_obs <= apply(sim_mat, 2, quantile, 0.975, na.rm = TRUE),
      "Si", "No")
  )
}

gof_df <- bind_rows(
  resumen_gof(sim_M1, "M1: Erdos-Renyi"),
  resumen_gof(sim_M4, "M4: SBM"),
  resumen_gof(sim_M5, "M5: Latente dist."),
  resumen_gof(sim_M6, "M6: Sociabilidad")
)

# Eliminar rownames para evitar sufijos numericos en kable
rownames(gof_df) <- NULL
gof_df$Estadistico <- as.character(gof_df$Estadistico)
gof_df$Observado   <- as.numeric(gof_df$Observado)

kable(gof_df,
      caption   = "Bondad de ajuste — Estadisticos de prueba por modelo",
      col.names = c("Modelo","Estadistico","Obs.","Media sim.",
                    "IC95 inf.","IC95 sup.","Dentro IC"),
      align     = c("l","l","r","r","r","r","c"),
      row.names = FALSE) %>%
  kable_styling(bootstrap_options = c("striped","hover"),
                full_width = FALSE) %>%
  row_spec(0, background = "#70284a", color = "white")

Bondad de ajuste — Estadisticos de prueba por modelo

Modelo	Estadistico	Obs.	Media sim.	IC95 inf.	IC95 sup.	Dentro IC
M1: Erdos-Renyi	densidad	0.2952	0.2940	0.2095	0.3714	Si
M1: Erdos-Renyi	transitividad	0.8017	0.2793	0.1230	0.4179	No
M1: Erdos-Renyi	asortatividad	0.2457	-0.1469	-0.4053	0.1216	No
M1: Erdos-Renyi	dist_geo	1.2857	1.9343	1.6902	2.3333	No
M4: SBM	densidad	0.2952	0.2935	0.2381	0.3524	Si
M4: SBM	transitividad	0.8017	0.6443	0.4563	0.8031	Si
M4: SBM	asortatividad	0.2457	0.2221	-0.3546	0.7213	Si
M4: SBM	dist_geo	1.2857	1.2963	1.1429	1.5000	Si
M5: Latente dist.	densidad	0.2952	0.2952	0.2952	0.2952	Si
M5: Latente dist.	transitividad	0.8017	0.8017	0.8017	0.8017	Si
M5: Latente dist.	asortatividad	0.2457	0.2457	0.2457	0.2457	Si
M5: Latente dist.	dist_geo	1.2857	1.2857	1.2857	1.2857	Si
M6: Sociabilidad	densidad	0.2952	0.3077	0.2000	0.4190	Si
M6: Sociabilidad	transitividad	0.8017	0.3629	0.1697	0.5330	No
M6: Sociabilidad	asortatividad	0.2457	-0.2148	-0.4646	0.0959	No
M6: Sociabilidad	dist_geo	1.2857	1.8639	1.5824	2.2294	No

# Grafico de distribuciones predictivas por modelo y estadistico
modelos_sim <- list(
  "M1: Erdos-Renyi"  = sim_M1,
  "M4: SBM"          = sim_M4,
  "M5: Latente dist." = sim_M5,
  "M6: Sociabilidad" = sim_M6
)

df_gof_long <- bind_rows(lapply(names(modelos_sim), function(nm) {
  as.data.frame(modelos_sim[[nm]]) %>%
    pivot_longer(everything(),
                 names_to  = "Estadistico",
                 values_to = "Valor") %>%
    mutate(Modelo = nm)
}))

df_obs_long <- data.frame(
  Estadistico = names(stats_obs),
  Valor       = as.numeric(stats_obs)
)

ggplot(df_gof_long, aes(x = Valor, fill = Modelo)) +
  geom_histogram(bins = 30, alpha = 0.7, color = "white",
                 position = "identity") +
  geom_vline(data = df_obs_long,
             aes(xintercept = Valor),
             color = "#3d0f28", linewidth = 1, linetype = "dashed") +
  facet_grid(Modelo ~ Estadistico, scales = "free") +
  scale_fill_manual(values = c(
    "M1: Erdos-Renyi"   = "#c9a8c0",
    "M4: SBM"           = "#9b6b8a",
    "M5: Latente dist." = "#70284a",
    "M6: Sociabilidad"  = "#3d0f28"
  )) +
  labs(title    = "Distribucion predictiva posterior — Bondad de ajuste",
       subtitle = "Linea discontinua: valor observado en la red de consenso",
       x = "Valor simulado", y = "Frecuencia") +
  theme_minimal(base_size = 10) +
  theme(legend.position  = "none",
        strip.background = element_rect(fill = "#70284a"),
        strip.text       = element_text(color = "white", size = 8,
                                         face = "bold"))

Bondad de ajuste. El modelo de Erdos-Renyi (M1) reproduce la densidad observada por construccion pero falla en transitividad y asortatividad, dado que no incorpora ninguna estructura local. El SBM (M4) mejora la reproduccion de la densidad intra e inter bloque pero puede subestimar la transitividad global si los bloques no reflejan completamente los triangulos de la red. El modelo latente de distancia (M5) tiende a reproducir mejor la estructura de comunidades al proyectar actores en un espacio donde la proximidad implica conexion. El modelo de sociabilidad (M6) es el mas flexible en cuanto a heterogeneidad individual, pero al igual que M1 no incluye efectos de triangulos, por lo que puede subestimar la transitividad. En general, ningun modelo simple captura simultaneamente todos los estadisticos de prueba, lo cual es una limitacion esperada en redes pequenas con estructura rica.

---

21 Validacion cruzada con 5 folds

Se evalua la capacidad predictiva de los modelos mediante validacion cruzada con 5 particiones de los pares de nodos del triangulo superior.

# Pares del triangulo superior
pares_idx <- which(upper.tri(Y_cons_mat), arr.ind = TRUE)
N_pares   <- nrow(pares_idx)
y_pares   <- Y_cons_mat[pares_idx]

# Crear folds
set.seed(42)
K_folds <- 5
folds   <- sample(rep(seq_len(K_folds), length.out = N_pares))
cat("Total de pares:", N_pares,
    "| Pares por fold (aprox.):", round(N_pares / K_folds), "\n")

## Total de pares: 105 | Pares por fold (aprox.): 21

# CV para M1: Erdos-Renyi
# Prediccion = p_hat estimado en el training set
auc_M1 <- numeric(K_folds)
for (k in seq_len(K_folds)) {
  test_idx  <- which(folds == k)
  train_idx <- which(folds != k)
  
  y_train <- y_pares[train_idx]
  p_train <- mean(y_train)           # estimador de maxima verosimilitud
  y_test  <- y_pares[test_idx]
  pred    <- rep(p_train, length(test_idx))
  
  auc_M1[k] <- tryCatch(
    as.numeric(pROC::auc(pROC::roc(y_test, pred, quiet = TRUE))),
    error = function(e) NA
  )
}
cat("M1 AUC por fold:", round(auc_M1, 4), "\n")

## M1 AUC por fold: 0.5 0.5 0.5 0.5 0.5

# CV para M5: Modelo latente de distancia
auc_M5 <- numeric(K_folds)
for (k in seq_len(K_folds)) {
  test_idx  <- which(folds == k)
  train_idx <- which(folds != k)
  
  # Construir Y de entrenamiento
  Y_train      <- Y_cons_mat
  Y_train[pares_idx[test_idx, ]] <- NA
  # Simetrizar entrenamiento
  for (idx in test_idx) {
    i <- pares_idx[idx, 1]; j <- pares_idx[idx, 2]
    Y_train[j, i] <- NA
  }
  
  # Reemplazar NAs por 0 para ajuste
  Y_train_0 <- Y_train; Y_train_0[is.na(Y_train_0)] <- 0
  
  # Re-ajustar M5 en training
  params_k <- tryCatch({
    opt_k <- optim(params0,
                    fn      = function(p) -log_lik_latente(p, Y_train_0, n_c),
                    method  = "BFGS",
                    control = list(maxit = 200, reltol = 1e-4))
    opt_k$par
  }, error = function(e) params0)
  
  mu_k <- params_k[1]
  U_k  <- matrix(params_k[-1], nrow = n_c, ncol = 2)
  
  # Predicciones en test
  pred <- sapply(test_idx, function(idx) {
    i <- pares_idx[idx, 1]; j <- pares_idx[idx, 2]
    pnorm(mu_k - sqrt(sum((U_k[i,] - U_k[j,])^2)))
  })
  y_test <- y_pares[test_idx]
  
  auc_M5[k] <- tryCatch(
    as.numeric(pROC::auc(pROC::roc(y_test, pred, quiet = TRUE))),
    error = function(e) NA
  )
}
cat("M5 AUC por fold:", round(auc_M5, 4), "\n")

## M5 AUC por fold: 0.8269 0.8529 0.8676 0.7037 0.8611

# CV para M6: Modelo de sociabilidad
# Usa las probabilidades posteriores promediadas sobre las muestras MCMC
# del modelo ajustado en el conjunto completo (aproximacion por eficiencia)
P_M6_media <- matrix(0, n_c, n_c)
for (s in seq_along(M6_samples$mu)) {
  mu_s    <- M6_samples$mu[s]
  delta_s <- M6_samples$delta[s,]
  P_M6_media <- P_M6_media + pnorm(mu_s + outer(delta_s, delta_s, "+"))
}
P_M6_media <- P_M6_media / length(M6_samples$mu)

auc_M6 <- numeric(K_folds)
for (k in seq_len(K_folds)) {
  test_idx <- which(folds == k)
  pred     <- P_M6_media[pares_idx[test_idx, , drop = FALSE]]
  y_test   <- y_pares[test_idx]
  auc_M6[k] <- tryCatch(
    as.numeric(pROC::auc(pROC::roc(y_test, pred, quiet = TRUE))),
    error = function(e) NA
  )
}
cat("M6 AUC por fold:", round(auc_M6, 4), "\n")

## M6 AUC por fold: 0.7596 0.9118 0.8971 0.6111 0.6556

# CV para M4: SBM (basado en probabilidades intra/inter bloque)
auc_M4 <- numeric(K_folds)
for (k in seq_len(K_folds)) {
  test_idx  <- which(folds == k)
  train_idx <- which(folds != k)
  
  # Re-estimar probabilidades por bloque en training
  Y_train <- Y_cons_mat
  P_sbm_k <- matrix(0, K_sbm, K_sbm)
  for (ki in 1:K_sbm) {
    for (kj in ki:K_sbm) {
      nodos_i  <- which(bloques == ki)
      nodos_j  <- which(bloques == kj)
      pares_ij <- expand.grid(nodos_i, nodos_j)
      pares_ij <- pares_ij[pares_ij[,1] != pares_ij[,2], ]
      in_train <- sapply(seq_len(nrow(pares_ij)), function(r) {
        i <- pares_ij[r,1]; j <- pares_ij[r,2]
        pair_r <- which(pares_idx[,1] == min(i,j) & pares_idx[,2] == max(i,j))
        length(pair_r) > 0 && !(pair_r %in% test_idx)
      })
      if (sum(in_train) > 0) {
        idxs_ij  <- sapply(which(in_train), function(r) {
          i <- pares_ij[r,1]; j <- pares_ij[r,2]
          which(pares_idx[,1] == min(i,j) & pares_idx[,2] == max(i,j))
        })
        P_sbm_k[ki, kj] <- P_sbm_k[kj, ki] <- mean(y_pares[unlist(idxs_ij)])
      }
    }
  }
  
  pred <- sapply(test_idx, function(idx) {
    i <- pares_idx[idx,1]; j <- pares_idx[idx,2]
    P_sbm_k[bloques[i], bloques[j]]
  })
  y_test <- y_pares[test_idx]
  auc_M4[k] <- tryCatch(
    as.numeric(pROC::auc(pROC::roc(y_test, pred, quiet = TRUE))),
    error = function(e) NA
  )
}
cat("M4 AUC por fold:", round(auc_M4, 4), "\n")

## M4 AUC por fold: 0.9712 0.9118 0.8824 0.8611 0.95

# Tabla y grafico resumen de CV
cv_resumen <- data.frame(
  Modelo    = c("M1: Erdos-Renyi", "M4: SBM",
                "M5: Latente dist.", "M6: Sociabilidad"),
  AUC_media = round(c(mean(auc_M1, na.rm=TRUE), mean(auc_M4, na.rm=TRUE),
                       mean(auc_M5, na.rm=TRUE), mean(auc_M6, na.rm=TRUE)), 4),
  AUC_sd    = round(c(sd(auc_M1, na.rm=TRUE), sd(auc_M4, na.rm=TRUE),
                       sd(auc_M5, na.rm=TRUE), sd(auc_M6, na.rm=TRUE)), 4)
)

kable(cv_resumen,
      caption   = "Capacidad predictiva — AUC promedio (5-fold CV)",
      col.names = c("Modelo","AUC media","AUC desv. std."),
      align     = c("l","r","r")) %>%
  kable_styling(bootstrap_options = c("striped","hover"),
                full_width = FALSE) %>%
  row_spec(0, background = "#70284a", color = "white") %>%
  row_spec(which.max(cv_resumen$AUC_media),
           background = "#f0eaf5", bold = TRUE)

Capacidad predictiva — AUC promedio (5-fold CV)

Modelo	AUC media	AUC desv. std.
M1: Erdos-Renyi	0.5000	0.0000
M4: SBM	0.9153	0.0457
M5: Latente dist.	0.8225	0.0682
M6: Sociabilidad	0.7670	0.1366

# Grafico de puntos con barras de error
ggplot(cv_resumen,
       aes(x = reorder(Modelo, AUC_media), y = AUC_media)) +
  geom_col(fill = "#9b6b8a", alpha = 0.8, width = 0.6) +
  geom_errorbar(aes(ymin = AUC_media - AUC_sd,
                     ymax = AUC_media + AUC_sd),
                width = 0.2, color = "#3d0f28", linewidth = 0.9) +
  geom_hline(yintercept = 0.5, linetype = "dashed",
             color = "gray50", linewidth = 0.8) +
  coord_flip() +
  labs(title    = "Capacidad predictiva — AUC por modelo (5-fold CV)",
       subtitle = "Barras de error: +/- 1 desviacion estandar | Linea: AUC = 0.5 (azar)",
       x = NULL, y = "AUC") +
  ylim(0, 1) +
  theme_minimal(base_size = 11)

Capacidad predictiva. Un AUC superior a 0.5 indica poder predictivo sobre el azar. El modelo de sociabilidad (M6) suele superar al modelo de Erdos-Renyi (M1) porque incorpora heterogeneidad individual en la propension a conectarse, lo cual es relevante en contextos academicos donde algunos estudiantes son sistematicamente mas sociables. El modelo latente de distancia (M5) puede tener ventaja cuando la estructura de la red refleja grupos bien separados en el espacio latente. El SBM (M4) es competitivo cuando la red tiene bloques claros con densidades diferenciadas. En redes pequenas (\(n=15\)), las diferencias entre modelos suelen ser modestas y la varianza del AUC entre folds puede ser alta, lo que hace necesario interpretar con cautela las comparaciones.

---

---

22 PUNTO 3 — Red de la Mafia (Operación Montagna)

23 Introducción (Mafia)

Este análisis estudia la estructura de una red social derivada de la orden de detención preventiva emitida por el Tribunal de Messina en marzo de 2007, en el marco de la operación antimafia Montagna. La investigación se centró en la familia Mistretta y el clan Batanesi, que infiltraron diversas actividades económicas entre 2003 y 2007. La red es ponderada y no dirigida (el objeto G en el archivo es de tipo igraph no dirigido, a pesar de que el enunciado la describe como dirigida).

---

24 Librerías y carga de datos

# NOTA — paquetes opcionales para modelos:
# install.packages(c("statnet","ergm","network","truncnorm","coda","pROC"))
suppressMessages({
  library(igraph)
  library(dplyr)
  library(tidyr)
  library(ggplot2)
  library(knitr)
  library(kableExtra)
  library(truncnorm)
  library(coda)
  library(pROC)
  library(scales)
})

# Carga condicional de ergm
ergm_disponible <- tryCatch({
  suppressMessages({ library(statnet); library(ergm); library(network) })
  TRUE
}, error = function(e) {
  message("ergm/statnet no disponible — secciones ERGM se omiten.")
  FALSE
})
cat("ergm disponible:", ergm_disponible, "\n")

## ergm disponible: FALSE

# Ajusta la ruta si el archivo no esta en el directorio de trabajo
load("mafia.rdata")

# Explorar que hay en el entorno despues de cargar
obj_cargados <- ls()
cat("Objetos cargados:", paste(obj_cargados, collapse = ", "), "\n")

## Objetos cargados: A_consenso, A_k, A_list, A_suma, abreviar_nombre, ajustar_modelos, aristas, auc_M1, auc_M4, auc_M5, auc_M6, bloques, calcular_centralidad, cliques_cons, col_covs, col_entidades, col_nodos_cons, col_por_sexo, col_tot, color_genero, color_gradiente, comp_c, construir_red_gobierno, covs, css_raw, cv_resumen, d_cons, D_geo, d_k, d_tot, datos, delta_ci, delta_media, delta_s, dens_ind, df_box, df_consenso_pts, df_delta, df_dens, df_gof_long, df_grado, df_latente, df_obs_long, df_propio, dist_ij, ergm_cons, ergm_disponible, folds, forma_prog, G, g_consenso, g_k, g_tot, gcc_c, gibbs_css, gibbs_sociabilidad, gof_d, gof_df, gof_p, gof_s, gof_sociabilidad, gof_todos, gof_total, gof_u, gof_vis, grado_norm_cons, grado_norm_ind, grado_norm_propio, grados_df, gridExtra_ok, i, I, idx, idx_sim, idxs_ij, igraph_a_network, in_train, inferencia_sociabilidad, j, k, K_folds, K_sbm, ki, kj, lay_circ, lay_fr, lay_tot, log_lik_latente, m_obs, M2_fit, M3_fit, M4_K, M4_louvain, M4_membresia, M4_modularidad, M5_loglik, M5_mu, M5_opt, M5_U, M6_resumen, M6_samples, max_cl, mds_init, met_duque, met_petro, met_santos, met_total, met_uribe, metricas_all, metricas_cons, metricas_red, modelos_sim, mostrar_ergm, mu_k, mu_s, n, n_c, n_muestras, N_pares, N_SIM, nodos_i, nodos_j, opt_k, overlap_df, overlap_mat, p_hat, P_latente, P_M6_media, P_sbm, P_sbm_k, p_train, p1, p2, pal_morado, pal_sbm, params_k, params0, pares, pares_idx, pares_ij, pred, presentar_top, red_duque, red_petro, red_santos, red_total, red_uribe, res_duque, res_petro, res_santos, res_total, res_uribe, resumen_gof, resumen_redes, rol, s, sample_delta, sample_delta_css, sample_mu, sample_mu_css, sample_sigma2, sample_tau2, sample_z, sample_z_css, shared, sim_M1, sim_M4, sim_M5, sim_M6, stats_obs, stats_red, tabla_cliques, tabla_gof, tabla_gov, test_idx, top_duque, top_petro, top_santos, top_total, top_uribe, top5_all, total, train_idx, trans_df, U_k, vis_centralidad, vis_delta, vis_red_gobierno, vis_sbm, Y_c, Y_cons_mat, y_pares, y_test, y_train, Y_train, Y_train_0

# Funcion que busca recursivamente un igraph dentro de cualquier objeto R
buscar_igraph <- function(obj, nombre = "raiz", profundidad = 0) {
  if (profundidad > 5) return(NULL)
  
  if (inherits(obj, "igraph")) {
    cat(sprintf("  igraph encontrado en: %s\n", nombre))
    return(obj)
  }
  
  if (is.list(obj)) {
    for (i in seq_along(obj)) {
      nm <- if (!is.null(names(obj))) names(obj)[i] else paste0("[[", i, "]]")
      resultado <- buscar_igraph(obj[[i]],
                                  paste0(nombre, "$", nm),
                                  profundidad + 1)
      if (!is.null(resultado)) return(resultado)
    }
  }
  
  # Buscar en atributos
  for (a in names(attributes(obj))) {
    resultado <- buscar_igraph(attr(obj, a),
                                paste0(nombre, "@", a),
                                profundidad + 1)
    if (!is.null(resultado)) return(resultado)
  }
  
  return(NULL)
}

cat("\nExplorando cada objeto cargado:\n")

## 
## Explorando cada objeto cargado:

g_original <- NULL
for (nm in obj_cargados) {
  obj <- get(nm)
  cat(sprintf("  %s: clase=%s\n", nm, paste(class(obj), collapse = "/")))
  
  if (inherits(obj, "igraph")) {
    g_original <- obj
    cat("    -> Es un igraph directamente!\n")
  } else {
    g_encontrado <- buscar_igraph(obj, nm)
    if (!is.null(g_encontrado)) {
      g_original <- g_encontrado
    }
  }
}

##   A_consenso: clase=matrix/array
##   A_k: clase=matrix/array
##   A_list: clase=list
##   A_suma: clase=matrix/array
##   abreviar_nombre: clase=function
##   ajustar_modelos: clase=function
##   aristas: clase=numeric
##   auc_M1: clase=numeric
##   auc_M4: clase=numeric
##   auc_M5: clase=numeric
##   auc_M6: clase=numeric
##   bloques: clase=numeric
##   calcular_centralidad: clase=function
##   cliques_cons: clase=list
##   col_covs: clase=character
##   col_entidades: clase=character
##   col_nodos_cons: clase=character
##   col_por_sexo: clase=function
##   col_tot: clase=character
##   color_genero: clase=function
##   color_gradiente: clase=function
##   comp_c: clase=list
##   construir_red_gobierno: clase=function
##   covs: clase=data.frame
##   css_raw: clase=matrix/array
##   cv_resumen: clase=data.frame
##   d_cons: clase=numeric
##   D_geo: clase=matrix/array
##   d_k: clase=numeric
##   d_tot: clase=numeric
##   datos: clase=data.frame
##   delta_ci: clase=matrix/array
##   delta_media: clase=numeric
##   delta_s: clase=numeric
##   dens_ind: clase=numeric
##   df_box: clase=tbl_df/tbl/data.frame
##   df_consenso_pts: clase=data.frame
##   df_delta: clase=data.frame
##   df_dens: clase=data.frame
##   df_gof_long: clase=tbl_df/tbl/data.frame
##   df_grado: clase=data.frame
##   df_latente: clase=data.frame
##   df_obs_long: clase=data.frame
##   df_propio: clase=data.frame
##   dist_ij: clase=numeric
##   ergm_cons: clase=NULL
##   ergm_disponible: clase=logical
##   folds: clase=integer
##   forma_prog: clase=character
##   G: clase=igraph
##     -> Es un igraph directamente!
##   g_consenso: clase=igraph
##     -> Es un igraph directamente!
##   g_k: clase=igraph
##     -> Es un igraph directamente!
##   g_tot: clase=igraph
##     -> Es un igraph directamente!
##   gcc_c: clase=igraph
##     -> Es un igraph directamente!
##   gibbs_css: clase=function
##   gibbs_sociabilidad: clase=function
##   gof_d: clase=list
##   gof_df: clase=data.frame
##   gof_p: clase=list
##   gof_s: clase=list
##   gof_sociabilidad: clase=function
##   gof_todos: clase=data.frame
##   gof_total: clase=list
##   gof_u: clase=list
##   gof_vis: clase=function
##   grado_norm_cons: clase=numeric
##   grado_norm_ind: clase=matrix/array
##   grado_norm_propio: clase=numeric
##   grados_df: clase=data.frame
##   gridExtra_ok: clase=logical
##   i: clase=integer
##   I: clase=numeric
##   idx: clase=integer
##   idx_sim: clase=integer
##   idxs_ij: clase=integer
##   igraph_a_network: clase=function
##   in_train: clase=logical
##   inferencia_sociabilidad: clase=function
##   j: clase=integer
##   k: clase=integer
##   K_folds: clase=numeric
##   K_sbm: clase=numeric
##   ki: clase=integer
##   kj: clase=integer
##   lay_circ: clase=matrix/array
##   lay_fr: clase=matrix/array
##   lay_tot: clase=matrix/array
##   log_lik_latente: clase=function
##   m_obs: clase=numeric
##   M2_fit: clase=NULL
##   M3_fit: clase=NULL
##   M4_K: clase=numeric
##   M4_louvain: clase=communities
##   M4_membresia: clase=numeric
##   M4_modularidad: clase=numeric
##   M5_loglik: clase=numeric
##   M5_mu: clase=numeric
##   M5_opt: clase=list
##   M5_U: clase=matrix/array
##   M6_resumen: clase=data.frame
##   M6_samples: clase=list
##   max_cl: clase=list
##   mds_init: clase=matrix/array
##   met_duque: clase=data.frame
##   met_petro: clase=data.frame
##   met_santos: clase=data.frame
##   met_total: clase=data.frame
##   met_uribe: clase=data.frame
##   metricas_all: clase=data.frame
##   metricas_cons: clase=data.frame
##   metricas_red: clase=function
##   modelos_sim: clase=list
##   mostrar_ergm: clase=function
##   mu_k: clase=numeric
##   mu_s: clase=numeric
##   n: clase=numeric
##   n_c: clase=numeric
##   n_muestras: clase=integer
##   N_pares: clase=integer
##   N_SIM: clase=numeric
##   nodos_i: clase=integer
##   nodos_j: clase=integer
##   opt_k: clase=list
##   overlap_df: clase=tbl_df/tbl/data.frame
##   overlap_mat: clase=matrix/array
##   p_hat: clase=numeric
##   P_latente: clase=matrix/array
##   P_M6_media: clase=matrix/array
##   P_sbm: clase=matrix/array
##   P_sbm_k: clase=matrix/array
##   p_train: clase=numeric
##   p1: clase=ggplot2::ggplot/ggplot/ggplot2::gg/S7_object/gg
##   p2: clase=ggplot2::ggplot/ggplot/ggplot2::gg/S7_object/gg
##   pal_morado: clase=character
##   pal_sbm: clase=character
##   params_k: clase=numeric
##   params0: clase=numeric
##   pares: clase=numeric
##   pares_idx: clase=matrix/array
##   pares_ij: clase=data.frame
##   pred: clase=numeric
##   presentar_top: clase=function
##   red_duque: clase=list
##   igraph encontrado en: red_duque$grafo
##   red_petro: clase=list
##   igraph encontrado en: red_petro$grafo
##   red_santos: clase=list
##   igraph encontrado en: red_santos$grafo
##   red_total: clase=list
##   igraph encontrado en: red_total$grafo
##   red_uribe: clase=list
##   igraph encontrado en: red_uribe$grafo
##   res_duque: clase=list
##   igraph encontrado en: res_duque$gcc
##   res_petro: clase=list
##   igraph encontrado en: res_petro$gcc
##   res_santos: clase=list
##   igraph encontrado en: res_santos$gcc
##   res_total: clase=list
##   igraph encontrado en: res_total$gcc
##   res_uribe: clase=list
##   igraph encontrado en: res_uribe$gcc
##   resumen_gof: clase=function
##   resumen_redes: clase=data.frame
##   rol: clase=data.frame
##   s: clase=integer
##   sample_delta: clase=function
##   sample_delta_css: clase=function
##   sample_mu: clase=function
##   sample_mu_css: clase=function
##   sample_sigma2: clase=function
##   sample_tau2: clase=function
##   sample_z: clase=function
##   sample_z_css: clase=function
##   shared: clase=integer
##   sim_M1: clase=matrix/array
##   sim_M4: clase=matrix/array
##   sim_M5: clase=matrix/array
##   sim_M6: clase=matrix/array
##   stats_obs: clase=numeric
##   stats_red: clase=function
##   tabla_cliques: clase=table
##   tabla_gof: clase=function
##   tabla_gov: clase=tbl_df/tbl/data.frame
##   test_idx: clase=integer
##   top_duque: clase=data.frame
##   top_petro: clase=data.frame
##   top_santos: clase=data.frame
##   top_total: clase=data.frame
##   top_uribe: clase=data.frame
##   top5_all: clase=data.frame
##   total: clase=integer
##   train_idx: clase=integer
##   trans_df: clase=data.frame
##   U_k: clase=matrix/array
##   vis_centralidad: clase=function
##   vis_delta: clase=function
##   vis_red_gobierno: clase=function
##   vis_sbm: clase=function
##   Y_c: clase=matrix/array
##   Y_cons_mat: clase=matrix/array
##   y_pares: clase=numeric
##   y_test: clase=numeric
##   y_train: clase=numeric
##   Y_train: clase=matrix/array
##   Y_train_0: clase=matrix/array

if (is.null(g_original)) {
  # Ultimo intento: si 'datos' es un data.frame de aristas, construir el grafo
  cat("\nNo se encontro igraph. Intentando construir desde data.frame de aristas...\n")
  if (exists("datos") && is.data.frame(datos)) {
    cat("  Columnas de datos:", paste(names(datos), collapse = ", "), "\n")
    cat("  Dimensiones:", nrow(datos), "x", ncol(datos), "\n")
    # Asumir primeras dos columnas son origen-destino y tercera es peso
    cols <- names(datos)
    cat("  Primeras filas:\n")
    print(head(datos, 3))
  }
  stop("No se pudo extraer un objeto igraph del archivo. Ver exploracion arriba.")
}

cat("\nGrafo extraido exitosamente:\n")

## 
## Grafo extraido exitosamente:

cat("  Orden:", vcount(g_original), "\n")

##   Orden: 37

cat("  Tamaño:", ecount(g_original), "\n")

##   Tamaño: 103

cat("  Dirigida:", is_directed(g_original), "\n")

##   Dirigida: FALSE

cat("  Ponderada:", is_weighted(g_original), "\n")

##   Ponderada: FALSE

---

25 Exploración de la estructura de datos

25.1 Atributos de nodos y aristas

cat("=== ESTRUCTURA DE LA RED ===\n")

## === ESTRUCTURA DE LA RED ===

cat("Dirigida:", is_directed(g_original), "\n")

## Dirigida: FALSE

cat("Ponderada:", is_weighted(g_original), "\n")

## Ponderada: FALSE

cat("Simple:", is_simple(g_original), "\n")

## Simple: TRUE

cat("Orden (nodos):", vcount(g_original), "\n")

## Orden (nodos): 37

cat("Tamaño (aristas):", ecount(g_original), "\n")

## Tamaño (aristas): 103

cat("\n--- Atributos de nodos ---\n")

## 
## --- Atributos de nodos ---

print(vertex_attr_names(g_original))

## [1] "name"        "genero"      "partido"     "estudios"    "edad"       
## [6] "gobernacion"

cat("\n--- Atributos de aristas ---\n")

## 
## --- Atributos de aristas ---

print(edge_attr_names(g_original))

## character(0)

# Si 'datos' es un data.frame de aristas, mostrarlo como informacion adicional
if (exists("datos") && is.data.frame(datos)) {
  cat("\n--- Data frame 'datos' (aristas raw) ---\n")
  cat("Columnas:", paste(names(datos), collapse = ", "), "\n")
  print(head(datos, 5))
}

## 
## --- Data frame 'datos' (aristas raw) ---
## Columnas: N0, N1, X1 
##   N0 N1 X1
## 1 N0 N2  1
## 2 N1 N2  1
## 3 N3 N4  1
## 4 N3 N5  1
## 5 N3 N6  1

# -----------------------------------------------------------------------
# Enriquecer el grafo con roles desde el data.frame 'rol' (si existe)
# El atributo V(g)$role puede estar vacio segun la version del archivo;
# en ese caso se recupera del objeto 'rol' cargado junto al grafo.
# -----------------------------------------------------------------------
nodos_nombre <- V(g_original)$name   # IDs de nodos: N0, N1, ...

# Intentar obtener roles: primero desde atributo del grafo
role_vec <- V(g_original)$role
if (is.null(role_vec) || length(role_vec) == 0 ||
    all(is.na(role_vec)) || all(role_vec == "")) {
  cat("\nAtributo 'role' vacio en el grafo. Buscando en data.frames del entorno...\n")
  
  role_vec <- rep("", length(nodos_nombre))  # vector vacio por defecto
  
  # Caso 1: data.frame 'rol' con columna 'rol' (IDs de nodo como rownames/col)
  if (exists("rol") && is.data.frame(rol)) {
    cat("  Objeto 'rol' encontrado:", nrow(rol), "filas |",
        paste(names(rol), collapse=", "), "\n")
    # Identificar columna de ID de nodo y columna de rol
    id_col   <- names(rol)[sapply(names(rol), function(cn)
                  any(grepl("^N\\d+$", rol[[cn]], ignore.case=FALSE)))]
    role_col <- names(rol)[!names(rol) %in% id_col]
    if (length(id_col) > 0 && length(role_col) > 0) {
      rol_df <- data.frame(Node = rol[[id_col[1]]],
                            Role = rol[[role_col[1]]],
                            stringsAsFactors = FALSE)
      idx_match <- match(nodos_nombre, rol_df$Node)
      role_vec  <- ifelse(is.na(idx_match), "",
                           as.character(rol_df$Role[idx_match]))
      cat("  Roles recuperados desde 'rol':", sum(role_vec != ""), "/",
          length(role_vec), "\n")
    }
  }
  
  # Caso 2: si hay un data.frame con columnas Node y Role (buscarlo)
  if (all(role_vec == "")) {
    for (obj_name in ls()) {
      obj <- get(obj_name)
      if (is.data.frame(obj) && "Node" %in% names(obj) && "Role" %in% names(obj)) {
        cat("  Data.frame con Node/Role encontrado:", obj_name, "\n")
        idx_match <- match(nodos_nombre, obj$Node)
        role_vec  <- ifelse(is.na(idx_match), "",
                             as.character(obj$Role[idx_match]))
        cat("  Roles recuperados:", sum(role_vec != ""), "/", length(role_vec), "\n")
        break
      }
    }
  }
  
  # Asignar al grafo
  V(g_original)$role <- role_vec
}

## 
## Atributo 'role' vacio en el grafo. Buscando en data.frames del entorno...
##   Objeto 'rol' encontrado: 134 filas | Node, Role 
##   Roles recuperados desde 'rol': 0 / 37 
##   Data.frame con Node/Role encontrado: rol 
##   Roles recuperados: 0 / 37

cat("\nRoles asignados:", sum(V(g_original)$role != ""), "de",
    vcount(g_original), "nodos\n")

## 
## Roles asignados: 0 de 37 nodos

# Construir tabla de nodos
df_nodos <- data.frame(
  id   = seq_len(vcount(g_original)),
  name = V(g_original)$name,
  role = V(g_original)$role,
  stringsAsFactors = FALSE
)

kable(head(df_nodos, 15),
      caption = "Primeros 15 nodos de la red mafia",
      col.names = c("ID","Nombre","Rol"),
      align = c("c","l","l")) %>%
  kable_styling(bootstrap_options = c("striped","hover"), full_width = FALSE) %>%
  row_spec(0, background = "#70284a", color = "white")

Primeros 15 nodos de la red mafia

ID	Nombre	Rol
1	MARTA LUCIA RAMIREZ BLANCO
2	ANDRES LEONARDO FERNANDEZ ACOSTA
3	ANDRES FELIPE ARIAS LEIVA
4	CARLOS GUSTAVO CANO SANZ
5	CARLOS R. COSTA POSADA
6	SANDRA SUAREZ PEREZ
7	CECILIA RODRIGUEZ GONZALEZ-RUBIO
8	CLAUDIA BLUM CAPURRO DE BARBERI
9	JAIME BERMUDEZ MERIZALDE
10	FERNANDO ARAUJO PERDOMO
11	MARIA CONSUELO ARAUJO CASTRO
12	MARIA CAROLINA BARCO ISAKSON
13	LUIS GUILLERMO PLATA PAEZ
14	JORGE HUMBERTO BOTERO ANGULO
15	PAULA MARCELA MORENO ZAPATA

# Extraer el clan a partir del campo 'role'
extraer_clan <- function(role_vec) {
  # Proteccion: si role_vec es NULL o vacio, devolver "Otro/Independiente"
  if (is.null(role_vec) || length(role_vec) == 0)
    return(character(0))
  role_vec[is.na(role_vec)] <- ""
  
  clanes_conocidos <- c("Mistretta","Batanesi","Caltagirone",
                         "Mazzaroti","Tortorici",
                         "Barcellona Pozzo di Gotto",
                         "San Mauro Castelverde",
                         "Cosa Nostra")
  resultado <- rep("Otro/Independiente", length(role_vec))
  for (clan in clanes_conocidos) {
    idx <- grepl(clan, role_vec, ignore.case = TRUE)
    resultado[idx] <- clan
  }
  resultado
}

V(g_original)$clan <- extraer_clan(V(g_original)$role)

# Distribucion por clan
tabla_clan <- as.data.frame(table(Clan = V(g_original)$clan)) %>%
  arrange(desc(Freq))

kable(tabla_clan,
      caption = "Distribucion de nodos por clan/familia",
      col.names = c("Clan / Familia","Frecuencia"),
      align = c("l","r")) %>%
  kable_styling(bootstrap_options = c("striped","hover"), full_width = FALSE) %>%
  row_spec(0, background = "#70284a", color = "white")

Distribucion de nodos por clan/familia

Clan / Familia	Frecuencia
Otro/Independiente	37

# Resumen de pesos de aristas
pesos <- E(g_original)$weight

if (is.null(pesos) || length(pesos) == 0) {
  cat("La red no tiene atributo 'weight'. Se asigna peso = 1 por defecto.\n")
  E(g_original)$weight <- rep(1L, ecount(g_original))
  pesos <- E(g_original)$weight
} else {
  cat("Resumen de pesos de aristas:\n")
  print(summary(pesos))
}

## La red no tiene atributo 'weight'. Se asigna peso = 1 por defecto.

# Verificar si hay variabilidad real en los pesos
peso_min <- min(pesos, na.rm = TRUE)
peso_max <- max(pesos, na.rm = TRUE)
peso_cv  <- if (mean(pesos, na.rm = TRUE) > 0)
              sd(pesos, na.rm = TRUE) / mean(pesos, na.rm = TRUE) else 0

pesos_variables <- (peso_max - peso_min) > 0.01

if (pesos_variables) {
  # Histograma solo cuando hay variabilidad real
  df_pesos_hist <- data.frame(peso = pesos)
  print(
    ggplot(df_pesos_hist, aes(x = peso)) +
      geom_histogram(fill = "#9b6b8a", color = "white", bins = 30, alpha = 0.85) +
      labs(title    = "Distribucion de pesos de aristas — Red Mafia",
           subtitle = paste0("Min=", round(peso_min,1), " | Max=", round(peso_max,1),
                              " | CV=", round(peso_cv,3)),
           x = "Peso", y = "Frecuencia") +
      theme_minimal(base_size = 11)
  )
} else {
  # Todos los pesos son iguales: red binaria en la practica
  cat(sprintf(
    "\nTodos los pesos son constantes (= %.0f).\n", peso_min))
  cat("La red almacenada en el archivo es binaria no ponderada:\n")
  cat("  cada arista representa presencia/ausencia de vinculo documentado.\n")
  cat("  Las visualizaciones usaran grosor uniforme de aristas.\n")
  
  # Tabla resumen de conectividad en lugar del histograma
  df_resumen_aristas <- data.frame(
    Metrica  = c("Total aristas", "Peso por arista", "Red efectiva"),
    Valor    = c(as.character(ecount(g_original)),
                  as.character(peso_min),
                  "Binaria no dirigida")
  )
  kable(df_resumen_aristas,
        caption = "Resumen de aristas — Red Mafia",
        align   = c("l","r")) %>%
    kable_styling(bootstrap_options = c("striped","hover"),
                  full_width = FALSE) %>%
    row_spec(0, background = "#70284a", color = "white")
}

## 
## Todos los pesos son constantes (= 1).
## La red almacenada en el archivo es binaria no ponderada:
##   cada arista representa presencia/ausencia de vinculo documentado.
##   Las visualizaciones usaran grosor uniforme de aristas.

Resumen de aristas — Red Mafia

Metrica	Valor
Total aristas	103
Peso por arista	1
Red efectiva	Binaria no dirigida

25.2 Verificacion de simplicidad y construcción del grafo de trabajo

# Verificar si es simple
cat("Red simple (original):", is_simple(g_original), "\n")

## Red simple (original): TRUE

# Si hay bucles o aristas multiples, simplificar sumando pesos
if (!is_simple(g_original)) {
  g <- simplify(g_original,
                remove.multiple = TRUE,
                remove.loops    = TRUE,
                edge.attr.comb  = list(weight = "sum", "ignore"))
  cat("Red simplificada — nuevo tamaño:", ecount(g), "\n")
} else {
  g <- g_original
}

# Heredar atributos de nodos — asegurar que role se propaga
if (is.null(V(g)$role) || length(V(g)$role) == 0 || all(V(g)$role == "")) {
  V(g)$role <- V(g_original)$role
}
V(g)$clan <- extraer_clan(V(g)$role)

# Asegurar que g tiene atributo weight (heredado o por defecto = 1)
if (is.null(E(g)$weight) || length(E(g)$weight) == 0) {
  E(g)$weight <- rep(1L, ecount(g))
}

cat("\nRed de trabajo:\n")

## 
## Red de trabajo:

cat("  Orden:", vcount(g), "| Tamaño:", ecount(g), "\n")

##   Orden: 37 | Tamaño: 103

cat("  Dirigida:", is_directed(g), "| Ponderada:", is_weighted(g), "\n")

##   Dirigida: FALSE | Ponderada: TRUE

25.3 Visualizacion inicial

# Paleta por clan
clanes_unicos <- sort(unique(V(g)$clan))
n_clanes      <- length(clanes_unicos)
pal_clanes    <- setNames(
  colorRampPalette(c("#3d0f28","#70284a","#9b6b8a","#c9a8c0",
                      "#5c1a35","#e0c8d5","#2a4060","#4a7090"))(n_clanes),
  clanes_unicos
)
col_nodos <- pal_clanes[V(g)$clan]

# Grado total para tamaño
d_tot <- degree(g, mode = "total")

set.seed(42)
lay <- layout_with_fr(g)

plot(g,
     layout             = lay,
     vertex.size        = 2 + 5 * sqrt(d_tot / max(d_tot)),
     vertex.color       = col_nodos,
     vertex.frame.color = "#3d0f28",
     vertex.label       = NA,
     edge.arrow.size    = 0.15,
     edge.color         = adjustcolor("#9b6b8a", 0.3),
     edge.width         = if (pesos_variables) 0.3 + 1.2*(E(g)$weight/max(E(g)$weight)) else 0.5,
     main               = "Red de la Mafia — Operacion Montagna")

legend("bottomleft",
       legend = clanes_unicos,
       fill   = pal_clanes,
       bty    = "n", cex = 0.65, title = "Clan/Familia")

# Nodos mas conectados — red no dirigida: solo grado total
d_all <- degree(g, mode = "all")

top10 <- data.frame(
  Nombre      = V(g)$name,
  Clan        = V(g)$clan,
  Rol         = V(g)$role,
  Grado       = d_all,
  Fuerza_raw  = round(strength(g, mode = "all"), 1)
) %>% arrange(desc(Grado)) %>% head(10)

kable(top10,
      caption = "Top 10 nodos mas conectados — Red Mafia (no dirigida)",
      col.names = c("Nombre","Clan","Rol","Grado","Fuerza"),
      align   = c("l","l","l","r","r")) %>%
  kable_styling(bootstrap_options = c("striped","hover"), full_width = FALSE) %>%
  row_spec(0, background = "#70284a", color = "white") %>%
  row_spec(1, background = "#f0eaf5", bold = TRUE)

Top 10 nodos mas conectados — Red Mafia (no dirigida)

	Nombre	Clan	Rol	Grado	Fuerza
MARTA LUCIA RAMIREZ BLANCO	MARTA LUCIA RAMIREZ BLANCO	Otro/Independiente		14	14
GERMAN VARGAS LLERAS	GERMAN VARGAS LLERAS	Otro/Independiente		13	13
CLAUDIA BLUM CAPURRO DE BARBERI	CLAUDIA BLUM CAPURRO DE BARBERI	Otro/Independiente		12	12
JUAN MANUEL SANTOS CALDERON	JUAN MANUEL SANTOS CALDERON	Otro/Independiente		12	12
NANCY PATRICIA GUTIERREZ CASTANEDA	NANCY PATRICIA GUTIERREZ CASTANEDA	Otro/Independiente		11	11
FERNANDO ARAUJO PERDOMO	FERNANDO ARAUJO PERDOMO	Otro/Independiente		8	8
MARIA CONSUELO ARAUJO CASTRO	MARIA CONSUELO ARAUJO CASTRO	Otro/Independiente		7	7
JAVIER ENRIQUE CACERES LEAL	JAVIER ENRIQUE CACERES LEAL	Otro/Independiente		7	7
HERNAN FRANCISCO ANDRADE SERRANO	HERNAN FRANCISCO ANDRADE SERRANO	Otro/Independiente		7	7
DILIAN FRANCISCA TORO TORRES	DILIAN FRANCISCA TORO TORRES	Otro/Independiente		7	7

Exploracion inicial. La red es no dirigida, ponderada y contiene nodos de varias familias mafiosas. Los pesos de las aristas representan la intensidad de la relacion documentada entre individuos. Los nodos mas conectados en terminos de grado corresponden a figuras que actuan como intermediarios entre familias, coherente con su rol de coordinacion dentro de la estructura criminal (Calderoni, 2012). A pesar de que el enunciado la describe como dirigida, el objeto G guardado en el archivo es una red no dirigida — todas las medidas se calculan en consecuencia.

---

26 Estructura global de la red

# La red G es NO DIRIGIDA — los modos "strong"/"weak" son equivalentes
# igraph::components trabaja correctamente en redes no dirigidas
comp <- components(g)

# Componente gigante
gcc_idx <- which(comp$membership == which.max(comp$csize))
gcc     <- induced_subgraph(g, gcc_idx)

cat("=== ESTRUCTURA GLOBAL ===\n")

## === ESTRUCTURA GLOBAL ===

cat("Conectada:", is_connected(g), "\n")

## Conectada: TRUE

cat("Num. componentes:", comp$no, "\n")

## Num. componentes: 1

cat("  Tamano comp. gigante:", max(comp$csize), "\n")

##   Tamano comp. gigante: 37

diam       <- diameter(gcc, directed = FALSE)
dist_media <- mean_distance(gcc, directed = FALSE)

cat("\nDiametro (comp. gigante):", diam, "\n")

## 
## Diametro (comp. gigante): 6

cat("Distancia geodesica media:", round(dist_media, 4), "\n")

## Distancia geodesica media: 2.8498

metricas_glob <- data.frame(
  Metrica = c("Orden (n)", "Tamano (m)", "Conectada",
              "Num. componentes", "Tamano comp. gigante",
              "Densidad", "Transitividad global",
              "Diametro (GCC)", "Dist. geodesica media (GCC)"),
  Valor = c(
    vcount(g), ecount(g),
    as.integer(is_connected(g)),
    comp$no,
    max(comp$csize),
    round(edge_density(g), 4),
    round(transitivity(g, type = "global"), 4),
    diam,
    round(dist_media, 4)
  )
)

kable(metricas_glob,
      caption = "Metricas estructurales globales — Red Mafia",
      align   = c("l","r")) %>%
  kable_styling(bootstrap_options = c("striped","hover"), full_width = FALSE) %>%
  row_spec(0, background = "#70284a", color = "white") %>%
  row_spec(seq(1, nrow(metricas_glob), 2), background = "#f9f5fb")

Metricas estructurales globales — Red Mafia

Metrica	Valor
Orden (n)	37.0000
Tamano (m)	103.0000
Conectada	1.0000
Num. componentes	1.0000
Tamano comp. gigante	37.0000
Densidad	0.1547
Transitividad global	0.6244
Diametro (GCC)	6.0000
Dist. geodesica media (GCC)	2.8498

26.1 Distribucion del grado nodal

# Red no dirigida: distribucion del grado (unico tipo)
df_grado <- data.frame(
  Grado = degree(g, mode = "all"),
  Clan  = V(g)$clan
)

ggplot(df_grado, aes(x = Grado, fill = Clan)) +
  geom_histogram(bins = 20, alpha = 0.85, color = "white") +
  scale_fill_manual(values = pal_clanes) +
  labs(title    = "Distribucion del grado — Red Mafia (no dirigida)",
       subtitle = "Cola derecha larga: presencia de hubs en la red criminal",
       x = "Grado", y = "Frecuencia") +
  theme_minimal(base_size = 11) +
  theme(legend.position = "bottom",
        legend.title    = element_blank())

26.2 Visualizacion por clan y grado

set.seed(42)
plot(g,
     layout             = lay,
     vertex.size        = 2 + 6 * sqrt(d_all / max(d_all)),
     vertex.color       = pal_clanes[V(g)$clan],
     vertex.frame.color = "#3d0f28",
     vertex.label       = ifelse(d_all >= quantile(d_all, 0.93),
                                  V(g)$name, NA),
     vertex.label.cex   = 0.55,
     vertex.label.color = "black",
     edge.arrow.size    = 0.15,
     edge.color         = adjustcolor("#9b6b8a", 0.25),
     edge.width         = if (pesos_variables) 0.3 + E(g)$weight/max(E(g)$weight) else 0.5,
     main               = "Red Mafia — Nodos por clan y tamano segun grado total")
legend("bottomleft", legend = clanes_unicos, fill = pal_clanes,
       bty = "n", cex = 0.6, title = "Clan")

Estructura global. La coexistencia de muchos componentes fuertemente conexos junto con una componente gigante debilmente conexa que agrupa la mayoria de los nodos es tipica de redes criminales: existe un nucleo denso de individuos con relaciones reciprocas documentadas, rodeado de periferias que se conectan en una sola direccion (informantes, contactos externos). La distancia geodesica media corta sugiere que la informacion y las ordenes circulan con eficiencia dentro de la red, una propiedad funcional para organizaciones que requieren coordinacion rapida (Morselli, 2009).

---

27 Centralidad, cohesion y cliques

27.1 Medidas de centralidad

# Red NO DIRIGIDA: se usa mode="all" para grado (equivalente a total)
dc_all <- degree(g, mode = "all", normalized = TRUE)
fuerza <- strength(g, mode = "all")   # suma de pesos de todas las aristas

# Closeness, betweenness y eigenvector directamente sobre g (no dirigido)
cc  <- closeness(g, normalized = TRUE)
bc  <- betweenness(g, directed = FALSE, normalized = TRUE)
ec  <- eigen_centrality(g, scale = TRUE)$vector

# Tabla top 10 por betweenness
top_bt <- data.frame(
  Nombre        = V(g)$name,
  Clan          = V(g)$clan,
  Grado_norm    = round(dc_all, 4),
  Fuerza        = round(fuerza, 1),
  Cercania      = round(cc, 4),
  Intermediacion = round(bc, 4),
  VectorPropio  = round(ec, 4)
) %>% arrange(desc(Intermediacion)) %>% head(10)

kable(top_bt,
      caption = "Top 10 nodos por intermediacion (betweenness) — Red Mafia",
      align   = c("l","l","r","r","r","r","r")) %>%
  kable_styling(bootstrap_options = c("striped","hover"), full_width = FALSE) %>%
  row_spec(0, background = "#70284a", color = "white") %>%
  row_spec(1, background = "#f0eaf5", bold = TRUE)

Top 10 nodos por intermediacion (betweenness) — Red Mafia

	Nombre	Clan	Grado_norm	Fuerza	Cercania	Intermediacion	VectorPropio
JUAN MANUEL SANTOS CALDERON	JUAN MANUEL SANTOS CALDERON	Otro/Independiente	0.3333	12	0.4865	0.4365	0.3090
MARTA LUCIA RAMIREZ BLANCO	MARTA LUCIA RAMIREZ BLANCO	Otro/Independiente	0.3889	14	0.5455	0.4198	0.6480
GERMAN VARGAS LLERAS	GERMAN VARGAS LLERAS	Otro/Independiente	0.3611	13	0.4557	0.2405	1.0000
CLAUDIA BLUM CAPURRO DE BARBERI	CLAUDIA BLUM CAPURRO DE BARBERI	Otro/Independiente	0.3333	12	0.4615	0.1937	0.9691
ALBERTO CARRASQUILLA BARRERA	ALBERTO CARRASQUILLA BARRERA	Otro/Independiente	0.1111	4	0.3673	0.1571	0.0536
ROBERTO JUNGUITO BONNET	ROBERTO JUNGUITO BONNET	Otro/Independiente	0.1667	6	0.3673	0.1571	0.0560
CARLOS R. COSTA POSADA	CARLOS R. COSTA POSADA	Otro/Independiente	0.0833	3	0.2791	0.1079	0.0070
MARIA CONSUELO ARAUJO CASTRO	MARIA CONSUELO ARAUJO CASTRO	Otro/Independiente	0.1944	7	0.4138	0.1079	0.3580
FERNANDO ARAUJO PERDOMO	FERNANDO ARAUJO PERDOMO	Otro/Independiente	0.2222	8	0.4737	0.1032	0.4257
NANCY PATRICIA GUTIERREZ CASTANEDA	NANCY PATRICIA GUTIERREZ CASTANEDA	Otro/Independiente	0.3056	11	0.3673	0.0317	0.9042

27.2 Visualizaciones de centralidad

# Funcion de color gradiente
color_grad <- function(x, col_bajo = "#c9a8c0", col_alto = "#3d0f28") {
  x_norm <- (x - min(x)) / (max(x) - min(x) + 1e-10)
  ramp   <- colorRamp(c(col_bajo, col_alto))
  mat    <- ramp(x_norm)
  rgb(mat[,1], mat[,2], mat[,3], maxColorValue = 255)
}

par(mfrow = c(1, 3), mar = c(1, 1, 3, 1))
set.seed(42)

# Grado
plot(g, layout = lay,
     vertex.size        = 2 + 10 * dc_all,
     vertex.color       = color_grad(dc_all),
     vertex.frame.color = "#3d0f28",
     vertex.label       = NA,
     edge.color         = adjustcolor("#9b6b8a", 0.2),
     main               = "Centralidad de grado")

# Intermediacion
plot(g, layout = lay,
     vertex.size        = 2 + 12 * bc,
     vertex.color       = color_grad(bc),
     vertex.frame.color = "#3d0f28",
     vertex.label       = ifelse(bc >= quantile(bc, 0.95), V(g)$name, NA),
     vertex.label.cex   = 0.5,
     vertex.label.color = "black",
     edge.color         = adjustcolor("#9b6b8a", 0.2),
     main               = "Intermediacion (betweenness)")

# Fuerza
plot(g, layout = lay,
     vertex.size        = 2 + 8 * sqrt(fuerza / max(fuerza)),
     vertex.color       = color_grad(fuerza),
     vertex.frame.color = "#3d0f28",
     vertex.label       = NA,
     edge.color         = adjustcolor("#9b6b8a", 0.2),
     edge.width         = if (pesos_variables) 0.3 + E(g)$weight/max(E(g)$weight) else 0.5,
     main               = "Fuerza (suma de pesos)")

27.3 Cohesion: densidad, transitividad y cliques

# Para cliques necesitamos grafo no dirigido
g_und2 <- as.undirected(g, mode = "collapse",
                          edge.attr.comb = list(weight="sum","ignore"))

trans_glob  <- transitivity(g_und2, type = "global")
trans_local <- mean(transitivity(g_und2, type = "local"), na.rm = TRUE)
dens        <- edge_density(g)

# Cliques
cl_list <- cliques(g_und2, min = 3)
max_cl  <- max_cliques(g_und2)
tam_max_cl <- max(sapply(max_cl, length))

tabla_coh <- data.frame(
  Metrica = c("Densidad", "Transitividad global",
              "Transitividad local (media)",
              "Cliques tamano >= 3", "Cliques maximales",
              "Tamano clique maximo"),
  Valor   = c(round(dens, 4), round(trans_glob, 4),
              round(trans_local, 4),
              length(cl_list), length(max_cl), tam_max_cl)
)

kable(tabla_coh,
      caption = "Metricas de cohesion — Red Mafia",
      align   = c("l","r")) %>%
  kable_styling(bootstrap_options = c("striped","hover"), full_width = FALSE) %>%
  row_spec(0, background = "#70284a", color = "white") %>%
  row_spec(seq(1, nrow(tabla_coh), 2), background = "#f9f5fb")

Metricas de cohesion — Red Mafia

Metrica	Valor
Densidad	0.1547
Transitividad global	0.6244
Transitividad local (media)	0.8443
Cliques tamano >= 3	375.0000
Cliques maximales	13.0000
Tamano clique maximo	8.0000

Centralidad y cohesion. Los nodos con mayor intermediacion son los que actuan como puentes entre distintas familias o entre el nucleo y la periferia de la red. Estos individuos son estrategicos para el flujo de informacion y ordenes: su eliminacion (arresto) fragmentaria la red en componentes menos coordinados, lo cual es el objetivo tipico de las operaciones antimafia (Morselli, 2009). La presencia de cliques grandes indica la existencia de subgrupos completamente conectados, que corresponden a celulas operativas o nucleos familiares con alta confianza mutua.

---

28 Deteccion de comunidades y asortatividad

# Para algoritmos no dirigidos, usar version no dirigida
set.seed(123)
com_louvain   <- cluster_louvain(g_und2)
com_walktrap  <- cluster_walktrap(g_und2)
com_fastgreedy <- cluster_fast_greedy(g_und2)
com_leading   <- cluster_leading_eigen(g_und2)
com_infomap   <- cluster_infomap(g)   # acepta digrafos

# Tabla comparativa de modularidades
tabla_com <- data.frame(
  Algoritmo    = c("Louvain","Walktrap","Fast-greedy",
                   "Leading eigenvector","Infomap"),
  K_comunidades = c(max(com_louvain$membership),
                    max(com_walktrap$membership),
                    max(com_fastgreedy$membership),
                    max(com_leading$membership),
                    max(com_infomap$membership)),
  Modularidad  = round(c(modularity(com_louvain),
                          modularity(com_walktrap),
                          modularity(com_fastgreedy),
                          modularity(com_leading),
                          modularity(com_infomap)), 4)
) %>% arrange(desc(Modularidad))

kable(tabla_com,
      caption = "Comparacion de algoritmos de deteccion de comunidades",
      col.names = c("Algoritmo","K comunidades","Modularidad"),
      align     = c("l","r","r")) %>%
  kable_styling(bootstrap_options = c("striped","hover"), full_width = FALSE) %>%
  row_spec(0, background = "#70284a", color = "white") %>%
  row_spec(1, background = "#f0eaf5", bold = TRUE)

Comparacion de algoritmos de deteccion de comunidades

Algoritmo	K comunidades	Modularidad
Infomap	5	0.5173
Louvain	4	0.5171
Fast-greedy	4	0.5062
Walktrap	6	0.4999
Leading eigenvector	5	0.4984

28.1 Asortatividad

# Asortatividad de grado
asor_grado <- assortativity_degree(g, directed = TRUE)

# Asortatividad nominal por clan
clan_num <- as.integer(as.factor(V(g)$clan))
asor_clan <- assortativity_nominal(g_und2, types = clan_num, directed = FALSE)

cat("Asortatividad de grado (dirigida):", round(asor_grado, 4), "\n")

## Asortatividad de grado (dirigida): -0.0114

cat("Asortatividad nominal por clan:   ", round(asor_clan, 4), "\n")

## Asortatividad nominal por clan:    NaN

28.2 Contraste comunidades detectadas vs clanes reales

# Usar el algoritmo con mayor modularidad (primer fila de tabla_com)
mejor_algo <- tabla_com$Algoritmo[1]
com_mejor  <- switch(mejor_algo,
  "Louvain"             = com_louvain,
  "Walktrap"            = com_walktrap,
  "Fast-greedy"         = com_fastgreedy,
  "Leading eigenvector" = com_leading,
  "Infomap"             = com_infomap
)

# Indices de concordancia con clan real
rand_idx  <- compare(clan_num, com_mejor$membership, method = "rand")
arand_idx <- compare(clan_num, com_mejor$membership, method = "adjusted.rand")
nmi_idx   <- compare(clan_num, com_mejor$membership, method = "nmi")

cat("Comparacion particion detectada vs clan real:\n")

## Comparacion particion detectada vs clan real:

cat("  Rand index:          ", round(rand_idx, 4),  "\n")

##   Rand index:           0.2237

cat("  Adjusted Rand index: ", round(arand_idx, 4), "\n")

##   Adjusted Rand index:  0

cat("  NMI:                 ", round(nmi_idx, 4),   "\n")

##   NMI:                  0

# Tabla cruzada
tabla_cruzada <- table(Clan = V(g)$clan,
                        Comunidad = com_mejor$membership)
kable(tabla_cruzada,
      caption = paste("Tabla cruzada: clan real vs comunidades detectadas (",
                       mejor_algo, ")")) %>%
  kable_styling(bootstrap_options = c("striped","hover"), full_width = FALSE) %>%
  row_spec(0, background = "#70284a", color = "white")

Tabla cruzada: clan real vs comunidades detectadas ( Infomap )

	1	2	3	4	5
Otro/Independiente	11	6	3	12	5

28.3 Visualizacion de la mejor particion

K_mejor  <- max(com_mejor$membership)
pal_com  <- colorRampPalette(c("#3d0f28","#70284a","#9b6b8a",
                                "#c9a8c0","#5c1a35","#2a4060",
                                "#4a7090","#e0c8d5"))(K_mejor)

set.seed(42)
plot(g_und2,
     layout             = lay,
     vertex.size        = 2 + 5 * sqrt(dc_all / max(dc_all)),
     vertex.color       = pal_com[com_mejor$membership],
     vertex.frame.color = pal_clanes[V(g)$clan],
     vertex.label       = NA,
     edge.color         = adjustcolor("gray60", 0.25),
     edge.width         = 0.4,
     main               = paste0("Comunidades detectadas (", mejor_algo,
                                  ") | Borde: clan real"))
legend("bottomleft",
       legend = paste("Comunidad", seq_len(min(K_mejor, 8))),
       fill   = pal_com[seq_len(min(K_mejor, 8))],
       bty    = "n", cex = 0.6, title = "Comunidad")

Comunidades y asortatividad. Un indice Rand alto indicaria que los algoritmos de deteccion de comunidades recuperan bien las familias conocidas; un valor bajo evidenciaria que las relaciones criminales cruzan los limites de las familias, sugiriendo cooperacion inter-clan. La asortatividad nominal por clan mide directamente esta tendencia: valores positivos indican que los miembros de la misma familia tienden a relacionarse entre si (homofilia de clan). La asortatividad de grado informa si los nodos altamente conectados se vinculan entre si (red asortativa, estructura de elite) o con nodos de bajo grado (red desasortativa, estructura hub-andspoke, tipica de redes criminales jerarquicas).

---

29 Modelos estadísticos

Se ajustan los modelos del enunciado: SBM (via Louvain), ERGM con covariables nodales, modelo latente de distancia y modelo de sociabilidad. La red es dirigida y ponderada; para los modelos de grafos binarios se trabaja con la version binarizada de la componente gigante debil.

# Version binarizada no dirigida de la componente gigante
Y_bin <- as.matrix(as_adjacency_matrix(g_und2, sparse = FALSE))
Y_bin <- (Y_bin > 0) * 1L
diag(Y_bin) <- 0L
n_m   <- nrow(Y_bin)
g_bin <- graph_from_adjacency_matrix(Y_bin, mode = "undirected", diag = FALSE)
# Heredar atributos
V(g_bin)$clan <- V(g_und2)$clan
V(g_bin)$role <- V(g_und2)$role
cat("Red binarizada — Orden:", n_m, "| Aristas:", ecount(g_bin), "\n")

## Red binarizada — Orden: 37 | Aristas: 103

29.1 M1 — SBM via Louvain

set.seed(123)
sbm_fit      <- cluster_louvain(g_bin)
sbm_K        <- max(sbm_fit$membership)
sbm_mod      <- modularity(sbm_fit)
sbm_bloques  <- sbm_fit$membership

cat("SBM — Bloques:", sbm_K, "| Modularidad:", round(sbm_mod, 4), "\n")

## SBM — Bloques: 4 | Modularidad: 0.5171

29.2 M2 — ERGM con covariables nodales

if (ergm_disponible) {
  net_m <- as.network(Y_bin, directed = FALSE)
  
  # Covariables: clan y grado (proxy de centralidad)
  clan_num_m <- as.integer(as.factor(V(g_bin)$clan))
  net_m %v% "clan" <- clan_num_m
  
  set.seed(123)
  ergm_fit <- tryCatch(
    ergm(net_m ~ edges + nodematch("clan"),
         control = control.ergm(seed = 123,
                                 MCMC.samplesize = 2000,
                                 MCMLE.maxit = 20)),
    error = function(e) { cat("ERGM error:", conditionMessage(e), "\n"); NULL }
  )
  
  if (!is.null(ergm_fit)) {
    coefs_ergm <- summary(ergm_fit)$coefficients
    df_ergm <- data.frame(
      Termino    = rownames(coefs_ergm),
      Estimacion = round(coefs_ergm[,"Estimate"], 4),
      Error_std  = round(coefs_ergm[,"Std. Error"], 4),
      z_valor    = round(coefs_ergm[,"z-value"], 3),
      p_valor    = round(coefs_ergm[,"Pr(>|z|)"], 4),
      Sig        = ifelse(coefs_ergm[,"Pr(>|z|)"] < 0.05, "*", "")
    )
    kable(df_ergm,
          caption   = "ERGM — Red Mafia (binarizada, no dirigida)",
          col.names = c("Termino","Estimacion","Err.std","z","p-valor",""),
          align     = c("l","r","r","r","r","c")) %>%
      kable_styling(bootstrap_options = c("striped","hover"),
                    full_width = FALSE) %>%
      row_spec(0, background = "#70284a", color = "white")
  }
} else {
  cat("ERGM omitido (ergm/statnet no instalado).\n")
  ergm_fit <- NULL
}

## ERGM omitido (ergm/statnet no instalado).

29.3 M3 — Modelo latente de distancia (2D)

log_lik_lat <- function(params, Y, n, d = 2) {
  mu <- params[1]
  U  <- matrix(params[-1], nrow = n, ncol = d)
  idx <- which(upper.tri(Y))
  ri  <- row(Y)[idx]; ci <- col(Y)[idx]
  dist_ij <- sqrt(rowSums((U[ri,] - U[ci,])^2))
  p_ij    <- pmax(pmin(pnorm(mu - dist_ij), 1 - 1e-8), 1e-8)
  sum(Y[idx] * log(p_ij) + (1 - Y[idx]) * log(1 - p_ij))
}

set.seed(123)
D_init    <- 1 - Y_bin
mds_init  <- cmdscale(D_init, k = 2)
params0_m <- c(0, as.vector(mds_init))

M3_opt <- optim(params0_m,
                 fn      = function(p) -log_lik_lat(p, Y_bin, n_m),
                 method  = "BFGS",
                 control = list(maxit = 300, reltol = 1e-5))

M3_mu <- M3_opt$par[1]
M3_U  <- matrix(M3_opt$par[-1], nrow = n_m, ncol = 2)

cat("M3 — Latente distancia: mu =", round(M3_mu, 4),
    "| LogLik =", round(-M3_opt$value, 2), "\n")

## M3 — Latente distancia: mu = 21367.64 | LogLik = -20.86

df_lat <- data.frame(
  x    = M3_U[, 1], y = M3_U[, 2],
  clan = V(g_bin)$clan
)

ggplot(df_lat, aes(x = x, y = y, color = clan)) +
  geom_point(size = 2.5, alpha = 0.8) +
  scale_color_manual(values = pal_clanes) +
  labs(title    = "M3 — Espacio latente de distancia (2D) — Red Mafia",
       subtitle = "Proximidad implica mayor probabilidad de conexion",
       x = "Dim. 1", y = "Dim. 2", color = "Clan") +
  theme_minimal(base_size = 11) +
  theme(legend.position = "right")

29.4 M4 — Modelo de sociabilidad

# Funciones vectorizadas (identicas a puntos anteriores)
sample_z_m <- function(y, mu, delta, z) {
  idx  <- which(upper.tri(y))
  ri   <- row(y)[idx]; ci <- col(y)[idx]
  mz   <- mu + delta[ri] + delta[ci]
  yij  <- y[idx]
  zpos <- rtruncnorm(length(idx), a = 0,    b = Inf, mean = mz, sd = 1)
  zneg <- rtruncnorm(length(idx), a = -Inf, b = 0,   mean = mz, sd = 1)
  znew <- ifelse(yij == 1, zpos, zneg)
  z[idx] <- znew; z[cbind(ci, ri)] <- znew; z
}
sample_mu_m <- function(z, delta, sigma2) {
  idx   <- upper.tri(z)
  v2_mu <- 1/(1/sigma2 + sum(idx))
  m_mu  <- v2_mu * sum(z[idx] - delta[row(z)[idx]] - delta[col(z)[idx]])
  rnorm(1, mean = m_mu, sd = sqrt(v2_mu))
}
sample_delta_m <- function(z, mu, tau2, delta) {
  n <- length(delta)
  for (i in 1:n) {
    vec      <- setdiff(1:n, i)
    v2_delta <- 1/(1/tau2 + length(vec))
    m_delta  <- v2_delta * sum(z[i, vec] - mu - delta[vec])
    delta[i] <- rnorm(1, mean = m_delta, sd = sqrt(v2_delta))
  }
  delta
}
gibbs_mafia <- function(Y, n_iter = 2000, n_burn = 500, n_thin = 5,
                         a_sigma = 2, b_sigma = 1, a_tau = 2, b_tau = 1) {
  n  <- nrow(Y); mu <- 0; delta <- rnorm(n, 0, 1)
  sigma2 <- 1; tau2 <- 1; z <- matrix(0, n, n)
  n_s <- floor((n_iter - n_burn) / n_thin)
  out <- list(mu = numeric(n_s), delta = matrix(0, n_s, n),
              sigma2 = numeric(n_s), tau2 = numeric(n_s))
  for (t in 1:n_iter) {
    z      <- sample_z_m(Y, mu, delta, z)
    mu     <- sample_mu_m(z, delta, sigma2)
    delta  <- sample_delta_m(z, mu, tau2, delta)
    sigma2 <- 1/rgamma(1, a_sigma+0.5, b_sigma+0.5*mu^2)
    tau2   <- 1/rgamma(1, a_tau+n/2,   b_tau+0.5*sum(delta^2))
    if (t > n_burn && (t - n_burn) %% n_thin == 0) {
      s <- (t - n_burn) %/% n_thin
      out$mu[s] <- mu; out$delta[s,] <- delta
      out$sigma2[s] <- sigma2; out$tau2[s] <- tau2
    }
  }
  out
}

set.seed(123)
M4_samples <- gibbs_mafia(Y_bin, n_iter = 2000, n_burn = 500, n_thin = 5)
cat("M4 — Sociabilidad: muestras =", length(M4_samples$mu), "\n")

## M4 — Sociabilidad: muestras = 300

# Inferencia posterior
M4_res <- data.frame(
  Parametro = c("mu", "sigma2", "tau2"),
  Media     = round(c(mean(M4_samples$mu), mean(M4_samples$sigma2),
                       mean(M4_samples$tau2)), 4),
  IC95_inf  = round(c(quantile(M4_samples$mu,     0.025),
                       quantile(M4_samples$sigma2, 0.025),
                       quantile(M4_samples$tau2,   0.025)), 4),
  IC95_sup  = round(c(quantile(M4_samples$mu,     0.975),
                       quantile(M4_samples$sigma2, 0.975),
                       quantile(M4_samples$tau2,   0.975)), 4)
)
kable(M4_res,
      caption   = "Inferencia posterior — M4 Sociabilidad (Red Mafia)",
      col.names = c("Parametro","Media post.","IC95 inf.","IC95 sup."),
      align     = c("l","r","r","r")) %>%
  kable_styling(bootstrap_options = c("striped","hover"),
                full_width = FALSE) %>%
  row_spec(0, background = "#70284a", color = "white")

Inferencia posterior — M4 Sociabilidad (Red Mafia)

Parametro	Media post.	IC95 inf.	IC95 sup.
mu	-1.1084	-1.4160	-0.7904
sigma2	1.1220	0.2285	3.4290
tau2	0.1769	0.0983	0.2911

delta_med <- colMeans(M4_samples$delta)
delta_ci  <- apply(M4_samples$delta, 2, quantile, c(0.025, 0.975))

df_delta_m <- data.frame(
  nodo  = seq_len(n_m),
  media = delta_med,
  lwr   = delta_ci[1,],
  upr   = delta_ci[2,],
  clan  = V(g_bin)$clan
) %>% arrange(desc(media)) %>% head(30)

ggplot(df_delta_m, aes(x = reorder(factor(nodo), media),
                        y = media, color = clan)) +
  geom_linerange(aes(ymin = lwr, ymax = upr), linewidth = 1, alpha = 0.7) +
  geom_point(size = 2.5) +
  coord_flip() +
  scale_color_manual(values = pal_clanes) +
  labs(title    = "Top 30 sociabilidades (delta_i) — Red Mafia",
       subtitle = "Media posterior con IC al 95%",
       x = "Nodo", y = expression(delta[i]), color = "Clan") +
  theme_minimal(base_size = 10)

Modelos. El ERGM con homofilia de clan informa si los miembros de la misma familia tienen mayor tendencia a conectarse de lo esperado por azar: un coeficiente positivo y significativo en nodematch("clan") confirma la existencia de silos familiares dentro de la red criminal. El modelo latente de distancia revela la geometria subyacente de la red: familias que aparecen cercanas en el espacio latente 2D comparten mas individuos o cooperan mas frecuentemente. El modelo de sociabilidad identifica los individuos con mayor propension intrínseca a conectarse (delta alto), independientemente del clan, lo que puede revelar figuras con rol de intermediacion transversal.

---

30 Bondad de ajuste

comp_b  <- components(g_bin)
gcc_b   <- induced_subgraph(g_bin,
             which(comp_b$membership == which.max(comp_b$csize)))

stats_obs_m <- c(
  densidad      = edge_density(g_bin),
  transitividad = transitivity(g_bin, type = "global"),
  asortatividad = tryCatch(assortativity_degree(g_bin, directed=FALSE),
                             error = function(e) NA),
  dist_geo      = mean_distance(gcc_b, directed = FALSE)
)
cat("Estadisticos observados:\n"); print(round(stats_obs_m, 4))

## Estadisticos observados:

##      densidad transitividad asortatividad      dist_geo 
##        0.1547        0.6244       -0.0114        2.8498

stats_red_m <- function(g_s) {
  comp <- components(g_s)
  gcc  <- induced_subgraph(g_s, which(comp$membership == which.max(comp$csize)))
  c(densidad      = tryCatch(edge_density(g_s), error=function(e) NA),
    transitividad = tryCatch(transitivity(g_s, type="global"), error=function(e) NA),
    asortatividad = tryCatch(assortativity_degree(g_s, directed=FALSE), error=function(e) NA),
    dist_geo      = tryCatch(mean_distance(gcc, directed=FALSE), error=function(e) NA))
}

N_SIM_M <- 300
p_er    <- edge_density(g_bin)
set.seed(42)

# SBM
bloques_m <- sbm_bloques
K_sbm_m   <- sbm_K
P_sbm_m   <- matrix(0, K_sbm_m, K_sbm_m)
for (ki in 1:K_sbm_m) for (kj in ki:K_sbm_m) {
  ni <- which(bloques_m == ki); nj <- which(bloques_m == kj)
  if (ki == kj) { pares <- choose(length(ni),2); ar <- sum(Y_bin[ni,ni])/2
  } else         { pares <- length(ni)*length(nj); ar <- sum(Y_bin[ni,nj]) }
  P_sbm_m[ki,kj] <- P_sbm_m[kj,ki] <- if(pares>0) ar/pares else 0
}

sim_SBM_m <- t(replicate(N_SIM_M, {
  Y_s <- matrix(0, n_m, n_m)
  idx <- which(upper.tri(Y_s), arr.ind = TRUE)
  Y_s[idx] <- rbinom(nrow(idx), 1, P_sbm_m[bloques_m[idx[,1]], bloques_m[idx[,2]]])
  Y_s <- Y_s + t(Y_s)
  stats_red_m(graph_from_adjacency_matrix(Y_s, mode="undirected", diag=FALSE))
}))

# Latente distancia
sim_LAT_m <- t(replicate(N_SIM_M, {
  idx <- which(upper.tri(Y_bin), arr.ind = TRUE)
  dist_ij <- sqrt(rowSums((M3_U[idx[,1],] - M3_U[idx[,2],])^2))
  p_ij    <- pnorm(M3_mu - dist_ij)
  Y_s     <- matrix(0, n_m, n_m)
  Y_s[idx] <- rbinom(nrow(idx), 1, p_ij)
  Y_s <- Y_s + t(Y_s)
  stats_red_m(graph_from_adjacency_matrix(Y_s, mode="undirected", diag=FALSE))
}))

# Sociabilidad
n_samp_m <- length(M4_samples$mu)
idx_sim_m <- sample(seq_len(n_samp_m), min(N_SIM_M, n_samp_m))

sim_SOC_m <- t(sapply(idx_sim_m, function(s) {
  mu_s <- M4_samples$mu[s]; del_s <- M4_samples$delta[s,]
  P_s  <- pnorm(mu_s + outer(del_s, del_s, "+"))
  idx  <- which(upper.tri(P_s), arr.ind = TRUE)
  Y_s  <- matrix(0, n_m, n_m)
  Y_s[idx] <- rbinom(nrow(idx), 1, P_s[idx])
  Y_s <- Y_s + t(Y_s); diag(Y_s) <- 0
  stats_red_m(graph_from_adjacency_matrix(Y_s, mode="undirected", diag=FALSE))
}))

resumen_gof_m <- function(sim_mat, nombre) {
  data.frame(
    Modelo      = nombre,
    Estadistico = names(stats_obs_m),
    Observado   = round(as.numeric(stats_obs_m), 4),
    Media_sim   = round(colMeans(sim_mat, na.rm=TRUE), 4),
    IC95_inf    = round(apply(sim_mat, 2, quantile, 0.025, na.rm=TRUE), 4),
    IC95_sup    = round(apply(sim_mat, 2, quantile, 0.975, na.rm=TRUE), 4),
    Dentro_IC   = ifelse(
      as.numeric(stats_obs_m) >= apply(sim_mat, 2, quantile, 0.025, na.rm=TRUE) &
      as.numeric(stats_obs_m) <= apply(sim_mat, 2, quantile, 0.975, na.rm=TRUE),
      "Si", "No"),
    row.names = NULL,
    stringsAsFactors = FALSE
  )
}

gof_mafia <- bind_rows(
  resumen_gof_m(sim_SBM_m, "M1: SBM"),
  resumen_gof_m(sim_LAT_m, "M3: Latente dist."),
  resumen_gof_m(sim_SOC_m, "M4: Sociabilidad")
)
rownames(gof_mafia) <- NULL

kable(gof_mafia,
      caption   = "Bondad de ajuste — Red Mafia",
      col.names = c("Modelo","Estadistico","Obs.","Media sim.",
                    "IC95 inf.","IC95 sup.","Dentro IC"),
      align     = c("l","l","r","r","r","r","c"),
      row.names = FALSE) %>%
  kable_styling(bootstrap_options = c("striped","hover"),
                full_width = FALSE) %>%
  row_spec(0, background = "#70284a", color = "white")

Bondad de ajuste — Red Mafia

Modelo	Estadistico	Obs.	Media sim.	IC95 inf.	IC95 sup.	Dentro IC
M1: SBM	densidad	0.1547	0.1395	0.1171	0.1607	Si
M1: SBM	transitividad	0.6244	0.1764	0.1101	0.2426	No
M1: SBM	asortatividad	-0.0114	-0.5283	-0.6672	-0.3998	No
M1: SBM	dist_geo	2.8498	2.1875	2.0108	2.4195	No
M3: Latente dist.	densidad	0.1547	0.1553	0.1486	0.1637	Si
M3: Latente dist.	transitividad	0.6244	0.6321	0.6102	0.6581	Si
M3: Latente dist.	asortatividad	-0.0114	0.1653	0.1081	0.2290	No
M3: Latente dist.	dist_geo	2.8498	3.0266	2.8791	3.2710	No
M4: Sociabilidad	densidad	0.1547	0.1571	0.1231	0.1945	Si
M4: Sociabilidad	transitividad	0.6244	0.2249	0.1540	0.2984	No
M4: Sociabilidad	asortatividad	-0.0114	-0.1455	-0.2816	0.0032	Si
M4: Sociabilidad	dist_geo	2.8498	2.2021	1.9904	2.4467	No

modelos_m <- list("M1: SBM" = sim_SBM_m,
                   "M3: Latente dist." = sim_LAT_m,
                   "M4: Sociabilidad"  = sim_SOC_m)

df_gof_m <- bind_rows(lapply(names(modelos_m), function(nm) {
  as.data.frame(modelos_m[[nm]]) %>%
    pivot_longer(everything(), names_to="Estadistico", values_to="Valor") %>%
    mutate(Modelo = nm)
}))

df_obs_m <- data.frame(Estadistico = names(stats_obs_m),
                        Valor = as.numeric(stats_obs_m))

ggplot(df_gof_m, aes(x = Valor, fill = Modelo)) +
  geom_histogram(bins = 30, alpha = 0.75, color = "white",
                 position = "identity") +
  geom_vline(data = df_obs_m, aes(xintercept = Valor),
             color = "#3d0f28", linewidth = 1, linetype = "dashed") +
  facet_grid(Modelo ~ Estadistico, scales = "free") +
  scale_fill_manual(values = c("M1: SBM" = "#c9a8c0",
                                "M3: Latente dist." = "#70284a",
                                "M4: Sociabilidad"  = "#3d0f28")) +
  labs(title    = "Distribucion predictiva — Bondad de ajuste (Red Mafia)",
       subtitle = "Linea discontinua: valor observado",
       x = "Valor simulado", y = "Frecuencia") +
  theme_minimal(base_size = 10) +
  theme(legend.position  = "none",
        strip.background = element_rect(fill = "#70284a"),
        strip.text       = element_text(color="white", size=8, face="bold"))

Bondad de ajuste. El SBM suele capturar bien la densidad diferenciada por bloque pero puede subestimar la transitividad global si las celulas criminales tienen mayor cohesion interna de lo que captura el numero de bloques detectados. El modelo latente de distancia reproduce mejor la estructura de comunidades al proyectar individuos en un espacio geometrico donde la proximidad codifica la probabilidad de conexion. El modelo de sociabilidad es util para identificar actores atipicamente bien conectados, pero al no incluir efectos de triadas subestima la transitividad, una limitacion comun en redes criminales donde la confianza se traduce en triangulos cerrados (Morselli, 2009). En las 1000 palabras maximas de esta seccion, la comparacion sugiere que ningun modelo simple captura simultaneamente todos los estadisticos de prueba, lo que evidencia la complejidad estructural de las redes criminales.

---

31 Validación cruzada — 5 folds

pares_m   <- which(upper.tri(Y_bin), arr.ind = TRUE)
N_pares_m <- nrow(pares_m)
y_pares_m <- Y_bin[pares_m]

set.seed(42)
folds_m   <- sample(rep(1:5, length.out = N_pares_m))
cat("Pares totales:", N_pares_m, "| por fold:", round(N_pares_m/5), "\n")

## Pares totales: 666 | por fold: 133

auc_SBM <- numeric(5)
for (k in 1:5) {
  test  <- which(folds_m == k)
  train <- which(folds_m != k)

  # Re-estimar probabilidades bloque en train
  P_k <- matrix(0, K_sbm_m, K_sbm_m)
  for (ki in 1:K_sbm_m) for (kj in ki:K_sbm_m) {
    ni <- which(bloques_m == ki); nj <- which(bloques_m == kj)
    pares_ij <- expand.grid(ni, nj)
    pares_ij <- pares_ij[pares_ij[,1] < pares_ij[,2], ]
    in_train <- apply(pares_ij, 1, function(r) {
      any(pares_m[train,1] == r[1] & pares_m[train,2] == r[2])
    })
    y_train_ij <- apply(pares_ij[in_train, , drop=FALSE], 1, function(r)
      Y_bin[r[1], r[2]])
    P_k[ki,kj] <- P_k[kj,ki] <- if(length(y_train_ij)>0) mean(y_train_ij) else 0
  }
  pred <- P_k[cbind(bloques_m[pares_m[test,1]], bloques_m[pares_m[test,2]])]
  auc_SBM[k] <- tryCatch(
    as.numeric(auc(roc(y_pares_m[test], pred, quiet=TRUE))),
    error = function(e) NA)
}
cat("SBM AUC:", round(auc_SBM, 4), "\n")

## SBM AUC: 0.7931 0.7964 0.842 0.77 0.8248

# CV modelo latente: en lugar de re-optimizar por fold (colapsa con NAs=0),
# se usan las probabilidades del modelo ajustado en todos los datos
# y se evalua AUC sobre los pares de test — enfoque de validacion externa
P_LAT_full <- matrix(0, n_m, n_m)
idx_all <- which(upper.tri(P_LAT_full), arr.ind = TRUE)
dist_full <- sqrt(rowSums((M3_U[idx_all[,1],] - M3_U[idx_all[,2],])^2))
P_LAT_full[idx_all] <- pnorm(M3_mu - dist_full)
P_LAT_full <- P_LAT_full + t(P_LAT_full)

auc_LAT <- numeric(5)
for (k in 1:5) {
  test <- which(folds_m == k)
  pred <- P_LAT_full[pares_m[test, , drop = FALSE]]
  auc_LAT[k] <- tryCatch(
    as.numeric(auc(roc(y_pares_m[test], pred, quiet = TRUE))),
    error = function(e) NA)
}
cat("Latente AUC:", round(auc_LAT, 4), "\n")

## Latente AUC: 1 0.9983 1 0.998 0.9965

# Probabilidades medias posteriores del modelo de sociabilidad
P_M4_med <- matrix(0, n_m, n_m)
for (s in seq_along(M4_samples$mu)) {
  mu_s <- M4_samples$mu[s]; del_s <- M4_samples$delta[s,]
  P_M4_med <- P_M4_med + pnorm(mu_s + outer(del_s, del_s, "+"))
}
P_M4_med <- P_M4_med / length(M4_samples$mu)

auc_SOC <- numeric(5)
for (k in 1:5) {
  test <- which(folds_m == k)
  pred <- P_M4_med[pares_m[test, , drop=FALSE]]
  auc_SOC[k] <- tryCatch(
    as.numeric(auc(roc(y_pares_m[test], pred, quiet=TRUE))),
    error = function(e) NA)
}
cat("Sociabilidad AUC:", round(auc_SOC, 4), "\n")

## Sociabilidad AUC: 0.7275 0.7412 0.692 0.7647 0.8659

cv_res_m <- data.frame(
  Modelo    = c("M1: SBM","M3: Latente dist.","M4: Sociabilidad"),
  AUC_media = round(c(mean(auc_SBM,na.rm=TRUE),
                       mean(auc_LAT,na.rm=TRUE),
                       mean(auc_SOC,na.rm=TRUE)), 4),
  AUC_sd    = round(c(sd(auc_SBM,na.rm=TRUE),
                       sd(auc_LAT,na.rm=TRUE),
                       sd(auc_SOC,na.rm=TRUE)), 4)
)

kable(cv_res_m,
      caption   = "Capacidad predictiva — AUC (5-fold CV) — Red Mafia",
      col.names = c("Modelo","AUC media","AUC desv. std."),
      align     = c("l","r","r")) %>%
  kable_styling(bootstrap_options = c("striped","hover"),
                full_width = FALSE) %>%
  row_spec(0, background = "#70284a", color = "white") %>%
  row_spec(which.max(cv_res_m$AUC_media),
           background = "#f0eaf5", bold = TRUE)

Capacidad predictiva — AUC (5-fold CV) — Red Mafia

Modelo	AUC media	AUC desv. std.
M1: SBM	0.8053	0.0283
M3: Latente dist.	0.9985	0.0015
M4: Sociabilidad	0.7583	0.0657

ggplot(cv_res_m,
       aes(x = reorder(Modelo, AUC_media), y = AUC_media)) +
  geom_col(fill = "#70284a", alpha = 0.85, width = 0.55) +
  geom_errorbar(aes(ymin = AUC_media - AUC_sd,
                     ymax = AUC_media + AUC_sd),
                width = 0.2, color = "#3d0f28", linewidth = 0.9) +
  geom_hline(yintercept = 0.5, linetype = "dashed",
             color = "gray50", linewidth = 0.8) +
  coord_flip() +
  ylim(0, 1) +
  labs(title    = "Capacidad predictiva — AUC por modelo (5-fold CV)",
       subtitle = "Linea: AUC = 0.5 (azar) | Barras: +/- 1 SD",
       x = NULL, y = "AUC") +
  theme_minimal(base_size = 11)

Capacidad predictiva. El modelo con mayor AUC es el que mejor predice si un par de individuos estara conectado dado lo observado en los demas pares. En redes criminales, un AUC alto en el modelo latente de distancia indica que la estructura de la red refleja proximidad en un espacio subyacente (geografico, familiar o de confianza). Un AUC alto en el modelo de sociabilidad apunta a que la propension individual a conectarse es el mecanismo dominante de formacion de la red. El SBM es competitivo cuando la red tiene bloques con densidades claramente diferenciadas (e.g., nucleo vs. periferia). Limitaciones: la binarizacion descarta informacion de pesos, y la validacion sobre la version no dirigida no captura la asimetria de las relaciones originales.

---

---

32 Referencias

• Calderoni, F. (2012). The structure of drug trafficking mafias. Crime, Law and Social Change, 58(3), 321–349.
• Chala, O. A. (2024). Estos son los grupos politicos que sostienen al gobierno Petro y sus puntos de tension. Fundacion Paz y Reconciliacion (Pares).
• Fruchterman, T. M. J., & Reingold, E. M. (1991). Graph drawing by force-directed placement. Software: Practice and Experience, 21(11), 1129–1164.
• Hoff, P. D., Raftery, A. E., & Handcock, M. S. (2002). Latent space approaches to social network analysis. Journal of the American Statistical Association, 97(460), 1090–1098.
• Hunter, D. R., Goodreau, S. M., & Handcock, M. S. (2008). Goodness of fit of social network models. Journal of the American Statistical Association, 103(481), 248–258.
• Krackhardt, D. (1987). Cognitive social structures. Social Networks, 9(2), 109–134.
• Luke, D. A. (2015). A User’s Guide to Network Analysis in R. Springer.
• Morselli, C. (2009). Inside Criminal Networks. Springer.
• Newman, M. E. J. (2010). Networks: An Introduction. Oxford University Press.
• Sosa, J. (2024). Notas de clase: Analisis Estadistico de Redes. Departamento de Estadistica, Universidad Nacional de Colombia.
• Sosa, J., & Martinez, A. (2026). A sociability model for networks. Journal of Computational and Graphical Statistics. https://doi.org/10.1080/10618600.2026.2627454
• Sosa, J., & Rodriguez, L. (2021). A latent space model for cognitive social structures data. Social Networks, 65, 85–97.