1 Introducción

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).

2 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: datos, ergm_disponible, G, rol
# 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
    }
  }
}
##   datos: clase=data.frame
##   ergm_disponible: clase=logical
##   G: clase=igraph
##     -> Es un igraph directamente!
##   rol: clase=data.frame
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: 143
cat("  Tamaño:", ecount(g_original), "\n")
##   Tamaño: 325
cat("  Dirigida:", is_directed(g_original), "\n")
##   Dirigida: FALSE
cat("  Ponderada:", is_weighted(g_original), "\n")
##   Ponderada: TRUE

3 Exploración de la estructura de datos

3.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: TRUE
cat("Simple:", is_simple(g_original), "\n")
## Simple: TRUE
cat("Orden (nodos):", vcount(g_original), "\n")
## Orden (nodos): 143
cat("Tamaño (aristas):", ecount(g_original), "\n")
## Tamaño (aristas): 325
cat("\n--- Atributos de nodos ---\n")
## 
## --- Atributos de nodos ---
print(vertex_attr_names(g_original))
## [1] "name" "role"
cat("\n--- Atributos de aristas ---\n")
## 
## --- Atributos de aristas ---
print(edge_attr_names(g_original))
## [1] "weight"
# 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
# Primeras filas 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 N0
2 N1
3 N3 boss family ‘’Caltagirone’’
4 N11 boss Cosa Nostra in Messina
5 N12 member family ‘’Mistretta’’
6 N4 enterpreneur
7 N5
8 N28
9 N25 executive family ‘’Mistretta’’
10 N6
11 N7
12 N8
13 N10 enterpreneur
14 N13 enterpreneur
15 N14 boss family ‘’Mazzaroti’’ 
# Extraer el clan a partir del campo 'role'
# Los roles siguen el patron: "tipo family 'CLAN'" o "tipo"
extraer_clan <- function(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 107
Batanesi 17
Mistretta 10
Mazzaroti 3
Tortorici 2
Barcellona Pozzo di Gotto 1
Caltagirone 1
Cosa Nostra 1
San Mauro Castelverde 1 
# Resumen de pesos de aristas
pesos <- E(g_original)$weight
cat("Resumen de pesos de aristas:\n")
## Resumen de pesos de aristas:
print(summary(pesos))
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   1.000   1.000   1.000   1.935   2.000  13.000
df_pesos_hist <- data.frame(peso = pesos)
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 = "Cada arista representa intensidad de relacion entre individuos",
       x = "Peso", y = "Frecuencia") +
  theme_minimal(base_size = 11)

3.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
V(g)$clan <- extraer_clan(V(g)$role)

cat("\nRed de trabajo:\n")
## 
## Red de trabajo:
cat("  Orden:", vcount(g), "| Tamaño:", ecount(g), "\n")
##   Orden: 143 | Tamaño: 325
cat("  Dirigida:", is_directed(g), "| Ponderada:", is_weighted(g), "\n")
##   Dirigida: FALSE | Ponderada: TRUE

3.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         = 0.3 + 1.2 * (E(g)$weight / max(E(g)$weight)),
     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
N18 N18 Mistretta executive family ‘’Mistretta’’ 40 86
N47 N47 Batanesi deputy boss family ‘’Batanesi’’ 26 100
N68 N68 Batanesi executive family ‘’Batanesi’’ 25 55
N27 N27 Batanesi executive family ‘’Batanesi’’ 19 48
N61 N61 Mistretta executive family ‘’Mistretta’’ 18 42
N22 N22 Otro/Independiente Pharmacist-member 17 43
N12 N12 Mistretta member family ‘’Mistretta’’ 16 31
N29 N29 Otro/Independiente enterpreneur 16 43
N11 N11 Cosa Nostra boss Cosa Nostra in Messina 15 32
N45 N45 Batanesi member family ‘’Batanesi’’ 14 40

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.

4 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: FALSE
cat("Num. componentes:", comp$no, "\n")
## Num. componentes: 5
cat("  Tamano comp. gigante:", max(comp$csize), "\n")
##   Tamano comp. gigante: 134
diam       <- diameter(gcc, directed = FALSE)
dist_media <- mean_distance(gcc, directed = FALSE)

cat("\nDiametro (comp. gigante):", diam, "\n")
## 
## Diametro (comp. gigante): 12
cat("Distancia geodesica media:", round(dist_media, 4), "\n")
## Distancia geodesica media: 4.1809
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) 143.0000
Tamano (m) 325.0000
Conectada 0.0000
Num. componentes 5.0000
Tamano comp. gigante 134.0000
Densidad 0.0320
Transitividad global 0.2785
Diametro (GCC) 12.0000
Dist. geodesica media (GCC) 4.1809

4.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())

4.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         = 0.3 + E(g)$weight / max(E(g)$weight),
     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).

5 Centralidad, cohesion y cliques

5.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
N18 N18 Mistretta 0.2817 86 0.4006 0.4654 0.1972
N47 N47 Batanesi 0.1831 100 0.3300 0.1992 1.0000
N68 N68 Batanesi 0.1761 55 0.3359 0.1944 0.4623
N61 N61 Mistretta 0.1268 42 0.3065 0.1559 0.0637
N27 N27 Batanesi 0.1338 48 0.3464 0.1190 0.5122
N43 N43 Otro/Independiente 0.0634 26 0.3500 0.1187 0.3556
N11 N11 Cosa Nostra 0.1056 32 0.3228 0.0911 0.0538
N75 N75 Mistretta 0.0775 22 0.2891 0.0848 0.0256
N12 N12 Mistretta 0.1127 31 0.3213 0.0785 0.0535
N76 N76 Otro/Independiente 0.0493 16 0.3260 0.0656 0.0321

5.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         = 0.3 + E(g)$weight / max(E(g)$weight),
     main               = "Fuerza (suma de pesos)")

5.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.0320
Transitividad global 0.2785
Transitividad local (media) 0.6760
Cliques tamano >= 3 585.0000
Cliques maximales 132.0000
Tamano clique maximo 7.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.

6 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
Fast-greedy 10 0.5717
Louvain 10 0.5699
Leading eigenvector 11 0.5650
Walktrap 11 0.5554
Infomap 17 0.5536

6.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.2326
cat("Asortatividad nominal por clan:   ", round(asor_clan, 4), "\n")
## Asortatividad nominal por clan:    0.1451

6.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.4612
cat("  Adjusted Rand index: ", round(arand_idx, 4), "\n")
##   Adjusted Rand index:  0.0221
cat("  NMI:                 ", round(nmi_idx, 4),   "\n")
##   NMI:                  0.2078
# 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 ( Fast-greedy )
1 2 3 4 5 6 7 8 9 10
Barcellona Pozzo di Gotto 0 1 0 0 0 0 0 0 0 0
Batanesi 1 14 0 0 2 0 0 0 0 0
Caltagirone 0 0 1 0 0 0 0 0 0 0
Cosa Nostra 0 0 1 0 0 0 0 0 0 0
Mazzaroti 0 0 1 0 2 0 0 0 0 0
Mistretta 3 0 4 2 0 1 0 0 0 0
Otro/Independiente 29 19 15 17 2 16 3 2 2 2
San Mauro Castelverde 0 0 1 0 0 0 0 0 0 0
Tortorici 0 1 0 0 1 0 0 0 0 0

6.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).

7 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: 143 | Aristas: 325

7.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: 9 | Modularidad: 0.5944

7.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).

7.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 = 0.3999 | LogLik = -791.29
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")

7.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 -2.2418 -2.4209 -2.0863
sigma2 2.2387 0.4709 8.3973
tau2 0.2352 0.1671 0.3264 
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.

8 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.0320        0.2785       -0.2326        3.1128
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.0320 0.0210 0.0195 0.0225 No
M1: SBM transitividad 0.2785 0.0205 0.0000 0.0328 No
M1: SBM asortatividad -0.2326 -0.8720 -0.8856 -0.8563 No
M1: SBM dist_geo 3.1128 2.1616 2.0847 2.3003 No
M3: Latente dist. densidad 0.0320 0.0324 0.0295 0.0354 Si
M3: Latente dist. transitividad 0.2785 0.1802 0.1435 0.2178 No
M3: Latente dist. asortatividad -0.2326 0.2527 0.1431 0.3799 No
M3: Latente dist. dist_geo 3.1128 4.0484 3.7971 4.3676 No
M4: Sociabilidad densidad 0.0320 0.0317 0.0274 0.0363 Si
M4: Sociabilidad transitividad 0.2785 0.1112 0.0853 0.1406 No
M4: Sociabilidad asortatividad -0.2326 -0.1980 -0.2682 -0.1268 Si
M4: Sociabilidad dist_geo 3.1128 2.8430 2.6928 3.0281 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.

9 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: 10153 | por fold: 2031
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.8466 0.8485 0.8492 0.8642 0.8817
# 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: 0.943 0.9603 0.9401 0.9575 0.9528
# 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.9166 0.8693 0.8703 0.8709 0.8839
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.8580 0.0150
M3: Latente dist. 0.9507 0.0089
M4: Sociabilidad 0.8822 0.0201 
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.

10 Referencias

  • Calderoni, F. (2012). The structure of drug trafficking mafias: the Ndrangheta and cocaine. Crime, Law and Social Change, 58(3), 321–349.
  • 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