library(readr)
library(readxl)
library(dplyr)
library(stringr)
library(purrr)
library(igraph)
library(ggraph)
library(tidygraph)
library(kableExtra)
library(tidyverse)
library(scales)
library(htmltools)La base de datos poder disponible en la página web del curso
corresponde a una red bipartita que relaciona funcionarios públicos con
las entidades públicas en las que fueron nombrados a lo largo de seis
gobiernos de Colombia, Uribe 1, Uribe 2, Santos 1, Santos 2, Duque y
Petro. Específicamente, estos datos se representan mediante una matriz
bipartita \(\mathbf{A}=[a_{i,j}]\) de
tamaño \(n\times m\), donde \(n\) representa el número de funcionarios y
\(m\) el número de entidades públicas.
En esta matriz, \(a_{i,j}=1\) si el
funcionario \(i\) fue nombrado en la
entidad \(j\), y \(a_{i,j}=0\) en caso contrario. La base
también incluye covariables individuales, como EDAD,
GENERO, PARTIDO_POLITICO,
NIVEL_DE_ESTUDIOS, UNIVERSIDAD,
CARRERA y GOBERNACION. Esta información fue
diseñada y recolectada junto con el estudiante Johan Andrés Hernández
Sarmiento como parte de su trabajo de grado. La base tiene carácter
estrictamente confidencial, por lo que ninguna persona cuenta con
autorización para compartirla, reproducirla o divulgarla.
El objetivo es llevar a cabo un análisis exhaustivo de las redes de poder en Colombia a través de distintos gobiernos. En el artículo Estos son los grupos políticos que sostienen al gobierno Petro y sus puntos de tensión, los autores desarrollan un ejercicio similar, aunque a una escala mucho menor y con recursos teóricos, empíricos y metodológicos más limitados desde el punto de vista del análisis de redes. A pesar de su simplicidad técnica, este texto constituye un buen punto de partida para comprender el propósito general de este caso de estudio.
Filtre la base de datos usando la covariable
GOBERNACION, de manera que se conserven todos aquellos
funcionarios que hicieron parte del Gobierno Petro, incluso si también
participaron en otros gobiernos. Sea \(\mathbf{A}_P\) la matriz bipartita
resultante de este filtro, de tamaño \(n_P\times m\), donde \(n_P\) corresponde al número de funcionarios
del Gobierno Petro y \(m\) corresponde
al número de entidades públicas. A partir de esta matriz, construya la
matriz de adyacencia ponderada \(\mathbf{A}_P\mathbf{A}_P^T\) de tamaño
\(n_P\times n_P\). Luego, construya la
matriz \(\mathbf{Y}_P\) a partir de
esta matriz de adyacencia ponderada, transformándola en una matriz
binaria y fijando en cero los elementos de la diagonal principal. Repita
este procedimiento para los gobiernos Duque, Santos y Uribe. En los
casos de Santos y Uribe, combine los dos periodos de cada gobierno en
una sola matriz. Al finalizar el proceso, debe contar con cuatro redes
binarias no dirigidas, una para cada gobierno.
De manera paralela y comparativa, construya visualizaciones decoradas para cada red, aprovechando al máximo la información nodal disponible. Presente estas visualizaciones de forma estratégica, de modo que faciliten la comparación entre gobiernos y permitan identificar patrones comunes, diferencias estructurales y posibles cambios en la configuración de las redes de poder.
Calcule medidas de centralidad para identificar cuáles son los funcionarios más influyentes dentro de cada red. A partir de estos resultados, construya visualizaciones decoradas que incorporen la información de centralidad de manera clara e informativa. Presente las visualizaciones de forma estratégica para comparar los resultados entre gobiernos. Finalmente, interprete los hallazgos de manera profunda, pero concisa, destacando sus posibles implicaciones políticas.
Analice cada una de las redes a nivel local y estructural mediante métricas de distancia, cohesión, conectividad y agrupamiento. Emplee un conjunto amplio de métricas que permita caracterizar de manera detallada las propiedades principales de cada red. Presente los resultados de forma estratégica, de modo que se facilite la comparación entre gobiernos y se identifiquen similitudes, diferencias y cambios relevantes en la estructura de las redes de poder. Finalmente, interprete los hallazgos de manera profunda, pero concisa, destacando sus posibles implicaciones políticas.
Ajuste, de manera independiente para cada red, un modelo de bloques estocásticos, un modelo de grafos aleatorios exponenciales y un modelo de socialidad con covariables. En el caso del modelo de grafos aleatorios exponenciales, incluya como términos del predictor lineal el número de aristas y las covariables con efectos homofílicos. Para el modelo de socialidad, incorpore también covariables con efectos homofílicos. Presente todos los resultados de forma estratégica, de modo que se facilite la comparación entre gobiernos y se identifiquen diferencias relevantes en la estructura de las redes de poder. Finalmente, interprete los resultados de manera profunda, pero concisa, destacando sus posibles implicaciones políticas.
Evalúe la bondad de ajuste de los modelos ajustados en el numeral anterior utilizando como estadísticos de prueba la densidad, la transitividad, la asortatividad y la distancia geodésica promedio. Presente los resultados en tablas y gráficos, interpretando detalladamente los hallazgos, identificando qué modelos capturan mejor las propiedades observadas en la red y limitando el análisis a un máximo de 1000 palabras.
Repita los numerales b., c. y d. para la red completa, sin filtrar por gobierno.
Solución Punto 1
En el siguiente bloque de código se cargan las bases de datos, se obtienen las matrices bipartitas \(\mathbf{A}_g\) para cada gobierno y sus respectivas proyecciones. El ersultado son cuatro matrices de adyacencia asociadas a grafos binarios no dirigidos.
# Base general
poder <- read_excel("Datos/poder.xlsx", sheet = "GABINETE_LIMPIO") %>%
select(-c(PREDECESORES, FLAG_MINISTERIO, FLAG_VIVO, X))
# Nivel de estudios
poder <- poder %>%
mutate(NIVEL_DE_ESTUDIOS =
case_when(NIVEL_DE_ESTUDIOS %in% c("MAESTRIA", "ESPECIALIZACION", "DOCTORADO", "POSGRADO") ~ "POSGRADO",
NIVEL_DE_ESTUDIOS %in% c("BACHILLERATO", "SIN BACHILLERATO", "TECNICO") ~ "SIN PREGRADO",
NIVEL_DE_ESTUDIOS == "NA" ~ "NA",
TRUE ~ "PREGRADO")
)
presidentes <- list(uribe = "URIBE-1|URIBE-2",
santos = "SANTOS-1|SANTOS-2",
duque = "DUQUE",
petro = "PETRO")
# Variables que no pertenecen a la matriz bipartita
vars_personales <- c("NOMBRE", "EDAD", "GENERO", "PARTIDO_POLITICO", "NIVEL_DE_ESTUDIOS", "UNIVERSIDAD", "CARRERA")
# Función para construir matriz bipartita
crear_bipartita <- function(datos) {
A <- datos %>%
select(-all_of(vars_personales)) %>%
mutate(across(everything(), as.numeric)) %>%
as.matrix()
rownames(A) <- datos$NOMBRE
A
}
# Bases por presidente
bases <- map(presidentes, ~ poder %>%
filter(str_detect(GOBERNACION, .x)) %>%
select(-GOBERNACION) )
bases <- map(bases, ~ .x %>%
mutate(EDAD = as.numeric(EDAD)))
# Matrices bipartitas
A <- map(bases, crear_bipartita)
# Matrices de adyacencia binarias
Y <- map(A, \(M) {
Y <- (M %*% t(M) > 0) * 1
diag(Y) <- 0
Y
})
list2env(setNames(Y, paste0("Y_", names(Y))), .GlobalEnv)## <environment: R_GlobalEnv>
## $uribe
## [1] 56 56
##
## $santos
## [1] 85 85
##
## $duque
## [1] 58 58
##
## $petro
## [1] 83 83
Los gobiernos con más funcionarios públicos nombrados fueron los de Santos, con 85 funcionarios. Mientras que durante los gobiernos de Uribe esta cantidad fue la menor, solamente con 56 funcionarios.
vars_nodales <- c("EDAD", "GENERO", "PARTIDO_POLITICO", "NIVEL_DE_ESTUDIOS", "UNIVERSIDAD", "CARRERA")
# Construir las redes
redes <- imap(Y, function(mat, nombre_gobierno){
g <- graph_from_adjacency_matrix(mat, mode = "undirected",diag = FALSE)
# Atributos nodales
atributos <- bases[[nombre_gobierno]] %>%
select(NOMBRE, all_of(vars_nodales))
atributos <- atributos[match(V(g)$name, atributos$NOMBRE), ]
for(v in vars_nodales){
g <- set_vertex_attr(g, name = v, value = atributos[[v]])
}
g
})Nota: Se hace una breve limpieza de la base de datos desde el archivo excel. Esta limpieza se aplicó a las variables nodales PARTIDO_POLITICO y NIVEL_DE_ESTUDIOS. En el primer caso, se limpiaron nombres de partidos ambiguos como ‘Partido alianza verde’ cuando ya hay otro llamado ‘Alianza verde’. En el segundo caso, esta variable se transforma para que sea el “Mayor nivel de estudios alcanzado”, de forma que si un funcionario tenía la entrada ‘ESPECIALIZACIÓN, MAESTRIA, DOCTORADO’ esta se toma como ‘DOCTORADO’.
Así mismo, de las variables ‘UNIVERSIDAD’ y ‘CARRERA’ se tomaron solo los primeros datos para cada individuo, con el fin de facilitar la comparación entre profesiones.
# Colores de los partidos
colores_partidos <- c("PARTIDO LIBERAL" = "red",
"PARTIDO CONSERVADOR" = "dodgerblue3",
"PARTIDO DE LA U" = "darkorange2",
"CAMBIO RADICAL" = "blue4",
"CENTRO DEMOCRATICO" = "deepskyblue",
"PARTIDO VERDE OXIGENO" = "green4",
"COLOMBIA RENACIENTE" = "limegreen",
"EN MARCHA" = "green",
"ALIANZA VERDE" = "darkgreen",
"PRIMERO COLOMBIA" = "black",
"CREEMOS PAIS" = "magenta",
"GENTE EN MOVIMIENTO" = "olivedrab",
"PACTO HISTORICO" = "blueviolet",
"COLOMBIA HUMANA" = "darkorchid1",
"PROMOVER COLOMBIA" = "pink",
"SALVACION NACIONAL" = "aquamarine",
"NUEVO LIBERALISMO" = "firebrick",
"PARTIDO COMUNISTA" = "gold3",
"POLO DEMOCRATICO ALTERNATIVO" = "yellow",
"SOY PORQUE SOMOS" = "orangered",
"NA" = "grey60")
graficar_red <- function(red, titulo, algoritmo) {
set.seed(123)
g_tbl <- as_tbl_graph(red)
layout_poder <- create_layout(g_tbl, layout = "igraph", algorithm = algoritmo)
layout_poder$educacion <- case_when(layout_poder$NIVEL_DE_ESTUDIOS == "NA" ~ 0,
layout_poder$NIVEL_DE_ESTUDIOS == "SIN PREGRADO" ~ 1,
layout_poder$NIVEL_DE_ESTUDIOS == "PREGRADO" ~ 2,
layout_poder$NIVEL_DE_ESTUDIOS == "POSGRADO" ~ 3)
ggraph(layout_poder) +
geom_edge_link(color = "grey70",
alpha = 0.5,
width = 1,
end_cap = circle(0, "mm") ) +
geom_node_point(aes(size = EDAD,
color = PARTIDO_POLITICO,
shape = GENERO)
) +
scale_size(range = c(1, 5), guide = "none") +
scale_color_manual(values = colores_partidos,
name = "Partido político") +
scale_shape_discrete(name = "Género") +
theme_void() +
theme(legend.position = "right") +
labs(title = titulo,
subtitle = "El tamaño de los vértices es proporcional a la edad")
}Como primer análisis de las redes de poder generadas, hablemos de la distribución de genero, el cual se representa mediante la forma del nodo en los gráficos. La proporción de Hombres y Mujeres en el gabinete estuvo más cerca de ser 50-50 durante el gobierno de Gustavo Petro con un 39.7% de mujeres. Por otro l ado, esta repartición fue más dispareja durante los gobiernos de Juan Manuel Santos con solo un 20% de mujeres.
En cuanto a los partidos políticos de los funcionarios vemos patrones muy interesantes. Durante los mandatos de Álvaro Uribe, primaban los funcionarios partidarios de ideologías que hoy se consideran de derecha, centro-derecha o centro; los más representativos de esta etapa son el Centro democrático, Cambio radical, Conservador y el Partido de la U.
El siguiente presidente trajo consigo cierto cambio en este paradigma: el partido conservador se clasifica como el más frecuente seguido del Liberal y del Partido de la U, sin que se note una clara dominación de alguno (lo que si se notaba en los periodos anteriores con el Centro Democrático). Vemos además, que se incorporan funcionarios de corrientes ideológicas más variadas, como del partido Alianza verde, Colombia Humana y el Pacto Histórico.
En el 2018 Iván Duque de cierto modo retomó una repartición similar a la de Álvaro Uribe. El Centro democrático domina drásticamente con 17 funcionarios de este partido, acompañados de partidos como Cambio radical, Partido conservador (estos dos con 5 funcionarios cada uno) y el partido de la U. Note que los partidos de centro-izquierda e izquierda que estuvieron presentes entre 2010 y 2018 ya no hacen parte de este gabinete.
Por último, con la llegada de Gustavo Petro naturalmente se observa una dominación de la izquierda y centro-izquierda en su gabinete. El Pacto Histórico cuenta con 18 funcionarios, seguido de Colombia Humana y el Partido Liberal, ambos con 7 funcionarios. Así mismo, el Partido comunista, el Polo democrático y el Nuevo liberalismo tienen su primera representación en el gabinete presidencial desde 2002.
datos_educacion <- map_dfr(names(bases), ~ bases[[.x]] %>%
mutate(GOBERNACION = .x) ) %>%
mutate(GOBERNACION = factor(GOBERNACION, levels = names(bases) )
) %>%
select(NOMBRE, NIVEL_DE_ESTUDIOS, GOBERNACION)
datos_educacion <- datos_educacion %>%
mutate(NIVEL_DE_ESTUDIOS = factor(NIVEL_DE_ESTUDIOS, levels = c("NA", "SIN PREGRADO", "PREGRADO", "POSGRADO"))
)
ggplot(datos_educacion,
aes(x = GOBERNACION,
fill = NIVEL_DE_ESTUDIOS)
) +
geom_bar(position = "fill") +
scale_y_continuous(labels = scales::percent) +
labs(x = NULL,
y = "% de funcionarios",
fill = "Nivel educativo",
title = "Distribución educativa por gobierno") +
theme_minimal()De forma complementaria, se presenta la distribución de los niveles educativos más altos alcanzados por los funcionarios de cada gobierno (se presenta de esta manera para no saturar la visualización de la red). Note que solo en el gobierno de Iván Duque todos los funcionarios tenían al menos un título Universitario, además que en este gobierno se tiene la mayor proporción de personas con posgrado. Durante el gobierno de Gustavo Petro la proporción de funcionarios con Posgrado disminuye considerablemente.
obtener_top <- function(g, medida, n = 3){
valores <- switch(medida,
"Closeness" = closeness(g, normalized = TRUE),
"Betweenness" = betweenness(g, normalized = TRUE),
"Eigenvector" = eigen_centrality(g)$vector)
top <- head(sort(valores, decreasing = TRUE), n)
paste(names(top), collapse = "<br>")
}
centralidades <- data.frame(Gobierno = names(redes),
Closeness = NA, Betweenness = NA, Eigenvector = NA)
for(i in seq_along(redes)){
g <- redes[[i]]
centralidades$Closeness[i] <- obtener_top(g, "Closeness")
centralidades$Betweenness[i] <- obtener_top(g, "Betweenness")
centralidades$Eigenvector[i] <- obtener_top(g, "Eigenvector")
}
kable(centralidades, align = "c", escape = FALSE, format = "html",
caption = "Top 3 de actores más centrales en las redes políticas colombianas") %>%
kable_styling(bootstrap_options = c("bordered", "striped", "hover"),
full_width = TRUE) %>%
column_spec(2:4, width = "25em")| Gobierno | Closeness | Betweenness | Eigenvector |
|---|---|---|---|
| uribe |
EVERTH BUSTAMANTE GARCIA DANIEL ANDRES GARCIA ARIZABALETA HERNAN MARTINEZ TORRES |
JUAN MANUEL SANTOS CALDERON MARTA LUCIA RAMIREZ BLANCO GERMAN VARGAS LLERAS |
GERMAN VARGAS LLERAS CLAUDIA BLUM CAPURRO DE BARBERI NANCY PATRICIA GUTIERREZ CASTANEDA |
| santos |
CLARA LUZ ROLDAN GONZALEZ ANDRES BOTERO PHILLIPSBOURNE JAIRO RAUL CLOPATOFSKY GHISAYS |
GERMAN VARGAS LLERAS MAURICIO CARDENAS SANTAMARIA FEDERICO ALONSO RENJIFO VELEZ |
GERMAN VARGAS LLERAS JUAN FERNANDO CRISTO BUSTOS LUIS FERNANDO VELASCO CHAVES |
| duque |
RODOLFO ENRIQUE ZEA NAVARRO ANDRES RAFAEL VALENCIA PINZON CARLOS EDUARDO CORREA ESCAF |
MARTA LUCIA RAMIREZ BLANCO JOSE MANUEL RESTREPO ABONDANO CARLOS HOLMES TRUJILLO GARCIA |
CLAUDIA BLUM CAPURRO DE BARBERI NANCY PATRICIA GUTIERREZ CASTANEDA ERNESTO JOSE MACIAS TOVAR |
| petro |
FRANCIA ELENA MARQUEZ MINA ANGELA YESENIA OLAYA REQUENE YANNAI KADAMANI FONRODONA |
LAURA CAMILA SARABIA TORRES OSCAR MAURICIO LIZCANO ARANGO ALEXANDER LOPEZ MAYA |
OSCAR MAURICIO LIZCANO ARANGO ARMANDO ALBERTO BENEDETTI VILLANEDA MIGUEL ANGEL PINTO HERNANDEZ |
Durante los gobiernos de Álvaro Uribe destacan Marta Lucia Ramirez (ministra de defensa), Juan Manuel Santos (ministro de defensa) y Germán Vargas Lleras (presidente del senado) por su centralidad de intermediación, lo que quiere decir que son importantes ya que se encuentran entre varios pares de nodos. Así mismo, Claudia Blum (embajadora ante la ONU) y Nancy Gutierrez (presidenta del senado) destacan por centralidad propia, en cuanto a que se rodean y relacionan con otros actores importantes. El senador Germán Vargas es el único que en esta etapa se considera de los más importantes bajo dos criterios.
En los 8 años siguientes, German Vargas Lleras se mantiene como funcionario importante en las redes de poder (bajo los criterios de intermediación y centralidad propia) debido a sus cargos de Ministro del Interior, Ministro de Vivienda y posterior Vicepresidente en el segundo mandato de Juan Manuel Santos. Otros funcionarios importantes debido a que se rodeaban de otros personajes centrales fueron Armando Benedetti y Luis Fernando Vlesco (presidentes del senado). Curiosamente, los tres funcionarios más importantes en términos de centralidad por cercanía se desempeñaron como directores de Coldeportes, el hoy Ministerio del Deporte.
Durante el gobierno de Iván Duque vuelven a destacar figuras que sobresalieron durante los gobiernos de Álvaro Uribe. Marta Lucia Ramirez vuelve, esta vez como vicepresidenta y canciller, así como Claudia Blum (canciller) y Nancy Gutierrez (ministra del interior), las últimas dos importantes en cuanto a la centralidad propia. La importancia de Jose Manuel Restrepo (ministro de comercio y de hacienda) y Carlos Holmes Trujillo (canciller y ministro de defensa) nos dice que ,junto a Marta Lucia, son funcionarios que se encuentran entre muchos caminos. Resaltan, por otro lado, Arturo Char (presidente del senado), Rodolfo Zea y Andres Valencia (ministros de agricultura), y Carlos Correa (ministro de ambiente).
Finalmente, en los últimos 4 años se escuchan más nombres nuevos. Francia Marquez (vicepresidenta y ministra de igualdad) es importante en cuanto a que está cerca de muchos otros funcionarios; al igual que Angela Olaya (ministra de ciencia) y Yannai Kadamani (ministra de cultura). La centralidad por intermediación resalta a Laura Camila Sanabria (jefa del gabinete presidencial, directora del DAPRE y canciller), Oscar Lizcano (director del DAPRE y ministro de las TIC), y Alexander Lopez (director del DNP y jefe del despacho presidencial). Por último, los funcionarios que durante este periodo se relacionarion con más personajes importantes fueron Oscar Lizcano, Armando Benedetti (ministro del interior y embajador) y Miguel Angel Pinto (presidente del senado).
library(visNetwork)
library(scales)
convertir_hex <- function(colores){
hex <- rgb(t(col2rgb(colores)), maxColorValue = 255)
names(hex) <- names(colores)
hex
}
graficar_red_interactiva <- function(g, titulo, top_n = 3,
centralidad = c("eigenvector", "closeness", "betweenness") ){
centralidad <- match.arg(centralidad)
colores_partidos_hex <- convertir_hex(colores_partidos)
# medidas de centralidad
cent <- switch(centralidad,
eigenvector = eigen_centrality(g)$vector,
closeness = closeness(g, normalized = TRUE),
betweenness = betweenness(g, normalized = TRUE))
top_actores <- if(top_n > 0){
names(sort(cent, decreasing = TRUE))[1:top_n]
} else {character(0)}
# nodos
nodes <- data.frame(
id = V(g)$name,
label = V(g)$name,
partido = trimws(as.character(V(g)$PARTIDO_POLITICO)),
genero = V(g)$GENERO,
edad = V(g)$EDAD,
carrera = V(g)$CARRERA,
nivel_edu = V(g)$NIVEL_DE_ESTUDIOS,
centralidad = cent,
stringsAsFactors = FALSE
)
nodes$central <- nodes$id %in% top_actores
nodes$label <- ifelse(nodes$central, nodes$id, "")
nodes$color <- lapply(seq_len(nrow(nodes)), function(i){
color_i <- colores_partidos_hex[nodes$partido[i]]
list(background = color_i, border = if(nodes$central[i]) "black" else color_i) })
nodes$value <- scales::rescale(nodes$edad, to = c(10, 40))
nodes$shape <- dplyr::case_when(nodes$genero == "F" ~ "dot",
nodes$genero == "M" ~ "triangle",
TRUE ~ "square")
# Borde si es central
nodes$borderWidth <- ifelse(nodes$central, 5, 1)
nodes$shadow <- nodes$central
# Tooltip
nodes$title <- paste0("<b>", nodes$id, "</b><br>",
"Partido: ", nodes$partido, "<br>",
"Nivel educativo: ", nodes$nivel_edu, "<br>",
"Carrera: ", nodes$carrera, "<br>",
"Edad: ", nodes$edad, "<br>",
toupper(substr(centralidad, 1, 1)),
substr(centralidad, 2, nchar(centralidad)),
": ",
round(nodes$centralidad, 3))
nodes$font <- lapply(
nodes$central, function(x){
if(x){list(size = 20, face = "bold", strokeWidth = 2)
} else {list(size = 0)}
})
# Aristas
edges <- igraph::as_data_frame(g, what = "edges")
names(edges)[1:2] <- c("from", "to")
edges$width <- 4
# Título de subtítulo dinámico
nombre_cent <- switch(centralidad,
eigenvector = "centralidad propia",
closeness = "centralidad por cercanía",
betweenness = "centralidad por intermediación")
titulo_html <- if(top_n > 0){
paste0(titulo,
"<br><span style='font-size:14px; font-weight:normal;'>",
"El tamaño de los nodos es proporcional a la edad. ",
"Se resaltan los ", top_n, " actores con mayor ", nombre_cent, ".</span>")
} else {
paste0(titulo,
"<br><span style='font-size:14px; font-weight:normal;'>",
"El tamaño de los nodos es proporcional a la edad. ",
".</span>")
}
# Graficar
visNetwork(nodes,
edges,
width = "100%",
height = "850px",
main = titulo_html
) %>%
visEdges(smooth = FALSE,
color = list(color = "rgba(180,180,180,0.5)")
) %>%
visOptions(
highlightNearest = list(enabled = TRUE, hover = TRUE),
selectedBy = list(variable = "partido", multiple = FALSE)
) %>%
visInteraction(dragNodes = TRUE,
dragView = TRUE,
zoomView = TRUE,
navigationButtons = TRUE) %>%
visPhysics(solver = "forceAtlas2Based",
stabilization = list(enabled = TRUE, iterations = 1000) ) %>%
visEvents(stabilizationIterationsDone ="function () {this.setOptions({physics:false});}" ) %>%
visLayout(randomSeed = 123)
}graficar_red_interactiva(redes$uribe, "Red de poder en Colombia - Gobiernos de Uribe", centralidad = "eigenvector")graficar_red_interactiva(redes$santos, "Red de poder en Colombia - Gobiernos de Santos", centralidad = "eigenvector")Se calculan el diámetro de la red y la distancia geodésica promedio, para cada uno de los gobiernos Adicionalmente se presenta la distribución de distancias geodésicas de las redes.
library(tibble)
distancia_redes <- data.frame(
Gobierno = names(redes),
Diametro = sapply(redes, diameter),
Distancia_promedio = sapply(
redes, function(g){mean_distance(g, directed = FALSE, unconnected = TRUE)} )
)
rownames(distancia_redes) <- NULL
distancia_redes <- distancia_redes %>%
column_to_rownames("Gobierno") %>%
t() %>%
as.data.frame()
distancia_redes["Diametro", ] <- sprintf("%.0f", as.numeric(distancia_redes["Diametro", ]))
distancia_redes["Distancia_promedio", ] <- sprintf("%.3f", as.numeric(distancia_redes["Distancia_promedio", ]))
kable(distancia_redes, align = "c", caption = "Medidas de distancia geodésica por gobierno") %>%
kable_styling(bootstrap_options = c("bordered", "striped", "hover"),
full_width = FALSE)| uribe | santos | duque | petro | |
|---|---|---|---|---|
| Diametro | 6 | 7 | 7 | 7 |
| Distancia_promedio | 2.814 | 2.873 | 2.849 | 2.890 |
Todas las redes de poder cuentan con un diámetro que podemos considerar ‘mediano’, dada la poca cantidad de nodos presente en cada red. Particularmente, en los gobiernos de Álvaro Uribe los funcionarios más alejados necesitan de 6 aristas para conectarse. Para los demás gobiernos, sus funcionarios necesitan una arista más para establecer una conexión.
En promedio, los funcionarios públicos se conectan mediante 2 o 3 intermediarios para todas las redes. Esto nos dice que en términos de distanciamiento, los distintos grafos no parecen tener diferencia alguna.
distancias_df <- lapply(
names(redes),
function(nombre){
g <- redes[[nombre]]
caminos <- distance_table(g)$res
data.frame(Gobierno = nombre,
Distancia = seq_along(caminos),
Frecuencia = prop.table(caminos))
} ) %>%
bind_rows()
ggplot(distancias_df, aes(x = Distancia,
y = Frecuencia) ) +
geom_col(fill = "grey60") +
facet_wrap(~ Gobierno) +
coord_cartesian(ylim = c(0, 0.8)) +
theme_minimal(base_size = 13) +
labs(title = "Distribución de distancias geodésicas",
subtitle = "Redes de poder de los gobiernos colombianos",
x = "Distancia geodésica",
y = "Frecuencia relativa")Curiosamente, la distancia geodésica más frecuente no es 1 (lo que sucede en varias redes sociales), sino 2 o 3. Además que las distribuciones no sean acampanadas nos indica variabilidad en el distanciamiento qeu tienen los actores dentro de la red.
En este apartado, se encuentran los ‘cliques’ más grandes de las redes junto con su frecuencia, además se llevan a cabo los censos diádicos. Del mismo modo, se calculan la densidad y transitividad.
tam_max_clique <- data.frame(
Gobierno = names(redes),
Clique_maximo = sapply(redes, clique_num),
Nodos = sapply(redes, vcount)
)
# proporción respecto al tamaño de la red
tam_max_clique$Proporcion <- round(tam_max_clique$Clique_maximo / tam_max_clique$Nodos, 2)
rownames(tam_max_clique) <- NULL
tabla_clique <- tam_max_clique %>%
select(Gobierno, Clique_maximo, Proporcion) %>%
tibble::column_to_rownames("Gobierno") %>%
t() %>%
as.data.frame()
rownames(tabla_clique) <- c("Tamaño clan máximo", "Proporción sobre total de nodos")
kable(tabla_clique, align = "c", caption = "Tamaño y proporción del clan máximo por gobiernos") %>%
kable_styling(bootstrap_options = c("bordered", "striped", "hover"), full_width = FALSE
)| uribe | santos | duque | petro | |
|---|---|---|---|---|
| Tamaño clan máximo | 8.00 | 10.00 | 6.0 | 10.00 |
| Proporción sobre total de nodos | 0.14 | 0.12 | 0.1 | 0.12 |
Que el tamaño del clan máximo represente una proporción pequeña del total de nodos de la red sugiere que los vínculos institucionales se encuentran fragmentados en pequeños grupos aislados, mostrando que no existe un núcleo amplio de funcionarios conectados entre sí a través de cargos comunes.
Dado que las conexiones representan experiencias compartidas en entidades públicas, un clan máximo de tamaño k indica la existencia de un grupo de k funcionarios cuyas trayectorias administrativas presentan coincidencias institucionales mutuas. Los gobiernos de Santos y Petro exhiben los clanes máximos más grandes (10 funcionarios), mientras que el gobierno Duque presenta el más pequeño (6 funcionarios), lo que sugiere una menor concentración de trayectorias compartidas dentro de su élite política.
cohesion <- data.frame(Gobierno = names(redes),
Densidad = round(sapply(redes, edge_density), 3),
Transitividad = round(sapply(redes, transitivity), 3) )
rownames(cohesion) <- NULL
cohesion <- cohesion %>%
tibble::column_to_rownames("Gobierno") %>%
t() %>%
as.data.frame()
kable(cohesion, align = "c") %>%
kable_styling(bootstrap_options = c("bordered", "striped", "hover"), full_width = FALSE)| uribe | santos | duque | petro | |
|---|---|---|---|---|
| Densidad | 0.075 | 0.078 | 0.047 | 0.051 |
| Transitividad | 0.629 | 0.639 | 0.629 | 0.652 |
# Transitividades locales
trans_df <- lapply(names(redes), function(nombre){
g <- redes[[nombre]]
data.frame(Gobierno = nombre,
Transitividad = transitivity( g, type = "local", isolates = "zero") ) }) %>%
bind_rows()
# Gráfico
ggplot(trans_df, aes(x = Transitividad)) +
geom_histogram(bins = 15,
fill = "grey60",
color = "white") +
facet_wrap(~ Gobierno) +
theme_minimal(base_size = 13) +
labs(title = "Distribución de la transitividad local",
x = "Transitividad local",y = "Frecuencia")La red de proyección más densa resultó ser la de los mandatos de Juan Manuel Santos, seguida de la de Álvaro Uribe, aunque ninguna está cerca de convertirse en un clan. La menos densa es la de los funcionarios de Iván Duque.
En términos de transitividad, esta se encuentra entre 0.62 y 0.66 para todas las redes de poder. El gobierno más transitivo fue el de Gustavo Petro, lo que quiere decir que en estos últimos cuatro años si el funcionario A compartió puesto con el funcionario B y este último con el C, el funcionario A es más propenso a trabajar en la misma institución que el C. Es natural que estos valores no sean tan altos ya que la transitividad no es del todo común en puestos de trabajo. Complementariamente, los gráficos de distribución de las transitividades locales muestran una alta dispersión para todas las redes.
Se estudia la conectividad por vértices y aristas, la componente gigante y sus características de conectividad.
## uribe santos duque petro
## FALSE FALSE FALSE FALSE
## uribe santos duque petro
## 0 0 0 0
## uribe santos duque petro
## 0 0 0 0
Como se vio en las visualizaciones, ninguna red está conectada.
## uribe santos duque petro
## 10 5 16 13
tam_componentes <- lapply(componentes, function(comp){table(sapply(comp, vcount)) })
tam_componentes## $uribe
##
## 1 2 3 37
## 2 4 3 1
##
## $santos
##
## 1 2 3 12 67
## 1 1 1 1 1
##
## $duque
##
## 1 2 3 4 26
## 4 7 2 2 1
##
## $petro
##
## 1 2 3 4 5 49
## 1 5 2 3 1 1
La red de co-pertenencia institucional de Iván Duque es la que cuenta con un mayor número de componentes (16), lo cual se debe en parte a la cantidad de nodos aislados con la que esta cuenta. La red con la menor cantidad de componentes es, como se esperaba, la de Juan Manuel Santos, con un valor de 5. El componente conectado más grande para cada gobierno tiene el tamaño de 37, 67, 26 y 49 nodos en orden cronológico respectivamente. Estos análisis nos hablan nuevamente de la alta conectividad que tiene la red del segundo presidente estudiado.
grafos_conectados <- lapply(redes, function(g){
comps <- decompose(g)
comps[[which.max(sapply(comps, vcount))]]
})
conectividad_df <- data.frame(
Gobierno = names(grafos_conectados),
Conectividad_nodal = sapply(grafos_conectados, vertex_connectivity),
Conectividad_aristas = sapply(grafos_conectados, edge_connectivity),
Puntos_articulacion = sapply(grafos_conectados, function(g) length(articulation_points(g))),
Proporcion_articulacion = sapply(grafos_conectados, function(g){length(articulation_points(g)) / vcount(g)})
)
kable(conectividad_df, digits = 3, align = "c",
caption = "Medidas de conectividad sobre la componente gigante, por carrera") %>%
kable_styling(bootstrap_options = c("bordered", "striped", "hover"),
full_width = FALSE)| Gobierno | Conectividad_nodal | Conectividad_aristas | Puntos_articulacion | Proporcion_articulacion | |
|---|---|---|---|---|---|
| uribe | uribe | 1 | 1 | 5 | 0.135 |
| santos | santos | 1 | 1 | 1 | 0.015 |
| duque | duque | 1 | 1 | 6 | 0.231 |
| petro | petro | 1 | 1 | 6 | 0.122 |
Ahora, con el fin de ir más allá en el análisis de la conectividad de las redes de proyección, se toman de base las componentes conexas de cada presidente y sobre estas se hacen los análisis. Todos las componentes gigantes muestran una conectividad nodal y de aristas de 1, lo que indica que basta con remover un nodo o una arista de las redes para que estas se desconecten.
La diferencia en la conectividad de las componentes está en la cantidad de puntos de articulación. Solo 1 de los funcionarios de los gobiernos Santos puede desconectar la red, mientras que para gobiernos como el de Duque y Petro esta cantidad es de 6 vértices.
obtener_comunidades <- function(g){
list(FastGreedy = cluster_fast_greedy(g),
LeadingEigen = cluster_leading_eigen(g),
Walktrap = cluster_walktrap(g),
Louvain = cluster_louvain(g),
LabelProp = cluster_label_prop(g),
Optimal = cluster_optimal(g),
Infomap = cluster_infomap(g)
)
}
# aplicar a todas las carreras
comunidades <- lapply(redes, obtener_comunidades)
set.seed(123)
modularidades <- lapply(names(comunidades), function(nombre){
algs <- comunidades[[nombre]]
g <- redes[[nombre]]
data.frame(Gobierno = nombre,
Algoritmo = names(algs),
Modularidad = sapply(algs, modularity),
Clusters = sapply(algs, length) ) }) %>%
bind_rows()
tabla_modularidad <- modularidades %>%
select(Gobierno, Algoritmo, Modularidad) %>%
pivot_wider(
names_from = Algoritmo,
values_from = Modularidad
)
rownames(tabla_modularidad) <- NULL
tabla_clusters <- modularidades %>%
select(Gobierno, Algoritmo, Clusters) %>%
pivot_wider(
names_from = Algoritmo,
values_from = Clusters
)
rownames(tabla_clusters) <- NULL
kable(tabla_modularidad, digits = 3, align = "c", caption = "Modularidad por algoritmo y gobierno") %>%
kable_styling(bootstrap_options = c("bordered", "striped", "hover"), full_width = FALSE
)| Gobierno | FastGreedy | LeadingEigen | Walktrap | Louvain | LabelProp | Optimal | Infomap |
|---|---|---|---|---|---|---|---|
| uribe | 0.591 | 0.599 | 0.583 | 0.606 | 0.579 | 0.606 | 0.601 |
| santos | 0.613 | 0.594 | 0.647 | 0.650 | 0.620 | 0.660 | 0.647 |
| duque | 0.700 | 0.700 | 0.704 | 0.714 | 0.696 | 0.714 | 0.708 |
| petro | 0.663 | 0.658 | 0.671 | 0.688 | 0.683 | 0.689 | 0.686 |
mejor_algoritmo <- modularidades %>%
group_by(Gobierno) %>%
slice_max(Modularidad, n = 1)
mejor_algoritmo## # A tibble: 4 × 4
## # Groups: Gobierno [4]
## Gobierno Algoritmo Modularidad Clusters
## <chr> <chr> <dbl> <int>
## 1 duque Optimal 0.714 19
## 2 petro Optimal 0.689 18
## 3 santos Optimal 0.660 9
## 4 uribe Optimal 0.606 13
En la búsqueda de comunidades para cada uno de los grafos de proyección funcionario-funcionario se implementan siete algoritmos de agrupamiento para cada uno de los gobiernos. En términos generales los que mejor se comportaron fueron los algoritmos Optimal y Louvain, habiendo empates en la modularidad obtenida por estos dos algoritmos para las redes de Uribe y Duque.
La mejor clasificación resulta en agrupamientos de 13, 9, 19 y 18 clusters para uribe, santos, duque y petro, respectivamente.
library(patchwork)
graficar_comunidades <- function(g, clustering, nombre, algoritmo){
g_tbl <- as_tbl_graph(g) %>%
mutate(comunidad = factor(clustering$membership),
grado = degree(g))
set.seed(123)
layout <- create_layout(g_tbl, layout = "igraph", algorithm = "nicely")
ggraph(layout) +
geom_edge_link(alpha = 0.5,
colour = "gray70",
width = 1,
end_cap = circle(0, "mm") ) +
geom_node_point(aes(color = comunidad,
size = EDAD,
shape = GENERO)
) +
scale_size(range = c(1, 5),guide = "none") +
scale_shape_discrete(name = "Género") +
guides(color = guide_legend(ncol = 2)
) +
theme_void() +
labs(title = paste("Agrupamiento -", nombre),
subtitle = paste("Algoritmo:", algoritmo, "| Modularidad =", round(modularity(clustering), 3) ),
color = "Comunidad"
)
}
plots_comunidades <- lapply(names(redes), function(nombre){
algoritmo <- mejor_algoritmo$Algoritmo[mejor_algoritmo$Gobierno == nombre]
clustering <- comunidades[[nombre]][[algoritmo]]
graficar_comunidades(g = redes[[nombre]],
clustering = clustering,
nombre = nombre,
algoritmo = algoritmo)
})
wrap_plots(plots_comunidades, ncol = 2) A continuación,se explora si las características nodales de los funcionarios repercuten en la presencia o no de enlaces entre ellos.
vars_nominales <- c( "GENERO", "PARTIDO_POLITICO", "NIVEL_DE_ESTUDIOS", "UNIVERSIDAD", "CARRERA")
vars_numericas <- c("EDAD")
tabla_asort <- data.frame(Red = names(redes))
# Variables categóricas
for(var in vars_nominales){
tabla_asort[[var]] <- sapply(redes, function(g){
assortativity_nominal(g,
types = as.integer(factor(vertex_attr(g, var))),
directed = FALSE)
})
}
# Variables numéricas
for(var in vars_numericas){
tabla_asort[[var]] <- sapply(redes, function(g){
# Eliminar nodos con NA en la variable
g2 <- delete_vertices(g, V(g)[is.na(vertex_attr(g, var))])
assortativity(g2,
values = vertex_attr(g2, var),
directed = FALSE)
})
}
kable(tabla_asort, digits = 3, align = "c",
caption = "Coeficientes de asortatividad por atributo") %>%
kable_styling(bootstrap_options = c("bordered") )| Red | GENERO | PARTIDO_POLITICO | NIVEL_DE_ESTUDIOS | UNIVERSIDAD | CARRERA | EDAD |
|---|---|---|---|---|---|---|
| uribe | 0.140 | 0.059 | -0.069 | 0.047 | 0.103 | -0.050 |
| santos | 0.136 | 0.039 | 0.076 | 0.064 | 0.125 | 0.222 |
| duque | -0.016 | 0.010 | 0.313 | 0.010 | 0.137 | 0.029 |
| petro | 0.215 | 0.051 | 0.098 | 0.001 | 0.047 | 0.042 |
En términos generales observamos valores de asortatividades positivas, lo que indica que los funcionarios públicos tienden a relacionarse con otros que tienen características similares. Analicemos esto a fondo.
De forma transversal, las variables partido político, universidad y carrera presentan asortatividades positivas. Esto deja ver que sin importar el gobernante al mando, los funcionarios tienden a ocupar los mismos cargos institucionales con aquellos otros que tengan estas características similares. Comportamiento a esperar tratándose de puestos que requieren (en teoría) un perfil de candidato.
Durante los gobiernos de Álvaro Uribe, los funcionarios se relacionaron con otros que tuviesen distinta edad y nivel de estudios. En el caso de Juan Manuel Santos, sus funcionarios son más propensos a relacionarse con aquellos que tienen edad similar. En la red de Duque hay una diferencia, y es que hombres y mujeres no tienden a ocupar los mismos cargos ya que la asortatividad para la variable genero fue negativa. Finalmente, para la red del gobierno Petro la universidad es una variable que no aporta casi nada en las conexiones entre nodos (asortatividad cercana a 0).
Ahora se calculan de manera independiente para cada red, un modelo de bloques estocásticos, un modelo de grafos aleatorios exponenciales y un modelo de sociabilidad con covariables.
Bajo este modelo, los vértices de una red están asociados a diferentes clases, y la propensión de formar aristas entre nodos está influenciada por la clase.
library(blockmodels)
sbm_modelos <- lapply(Y, function(mat) {
modelo <- BM_bernoulli(membership_type = "SBM_sym",
adj = mat,
verbosity = 0,
plotting = "")
modelo$estimate()
return(modelo)
})resultados_sbm <- lapply(names(sbm_modelos), function(nombre){
mod <- sbm_modelos[[nombre]]
# Número óptimo de clusters
Q <- which.max(mod$ICL)
# Probabilidades de pertenencia
Z <- mod$memberships[[Q]]$Z
# Asignaciones a los clusters
labs <- apply(Z, 1, which.max)
# Tamaño relativo de comunidades
alpha <- prop.table(table(labs))
# Probabilidades de interacción
Pi <- mod$model_parameters[[Q]]$pi
list(nombre = nombre,
ICL = mod$ICL,
Q = Q,
Z = Z,
labs = labs,
alpha = alpha,
Pi = Pi)
})
names(resultados_sbm) <- names(sbm_modelos)## uribe santos duque petro
## 3 5 2 3
De acuerdo con la verosimilitud de clasificación integrada (ICL) el número de grupos/bloques óptimos para las redes de poder son 3, 5, 2 y 3, para ÁLvaro Uribe, Juan Manuel Santos, Iván Duque y Gustavo Petro, respectivamente. Note que para todas las redes, el modelo clasifica en un número pequeño de clusters.
par(mfrow = c(2,2))
for(nombre in names(resultados_sbm)){
ICL <- resultados_sbm[[nombre]]$ICL
Q <- resultados_sbm[[nombre]]$Q
plot(ICL, type = "b", pch = 16, xlab = "Q", ylab = "ICL", main = nombre)
abline(v = Q, col = "red", lty = 2)
}Ahora, se obtienen las probabilidades de pertenencia a las comunidades \(E(\mathbf{z_i|\mathbf{Y=y}})\). Estas nos van a permitir más adelante determinar las etiquetas de asignación.
## $uribe
## [,1] [,2] [,3]
## [1,] 0.001782531 0.996434938 0.001782531
## [2,] 0.001782531 0.001782531 0.996434938
## [3,] 0.001782531 0.001782531 0.996434938
##
## $santos
## [,1] [,2] [,3] [,4] [,5]
## [1,] 0.001172333 0.001172333 0.001172333 0.995310668 0.001172333
## [2,] 0.001172333 0.001172333 0.001172333 0.995310668 0.001172333
## [3,] 0.001172333 0.001172333 0.001172333 0.001172333 0.995310668
##
## $duque
## [,1] [,2]
## [1,] 0.998275862 0.001724138
## [2,] 0.001724138 0.998275862
## [3,] 0.998275862 0.001724138
##
## $petro
## [,1] [,2] [,3]
## [1,] 0.001447554 0.997349077 0.001203369
## [2,] 0.001203369 0.997593261 0.001203369
## [3,] 0.997593261 0.001203369 0.001203369
## $uribe
##
## 1 2 3
## 12 12 32
##
## $santos
##
## 1 2 3 4 5
## 12 28 21 13 11
##
## $duque
##
## 1 2
## 40 18
##
## $petro
##
## 1 2 3
## 26 47 10
# Tamaño relativo de comunidades
lapply(resultados_sbm, \(x)
round(sort(x$alpha, decreasing = TRUE),3)
)## $uribe
## labs
## 3 1 2
## 0.571 0.214 0.214
##
## $santos
## labs
## 2 3 4 1 5
## 0.329 0.247 0.153 0.141 0.129
##
## $duque
## labs
## 1 2
## 0.69 0.31
##
## $petro
## labs
## 2 1 3
## 0.566 0.313 0.120
Para la red de poder de los gobiernos Uribe hay un cluster en el que se agrupan el 57,1% de los funcionarios, mientras que el restante 42,9% de nodos se reparte por igual en los otros dos clusters. En cuanto al presidente Santos, su cluster más grande contiene aproximadamente al 32,9% de los funcionarios, mientras que el menor se comprende de un 12,9%. Para la red de Iván Duque, uno de sus dos clusters dobla en tamaño al otro, con un 69% de los nodos. Por último, 56,6% de los funcioanrios del gobierno Petro se agrupan juntos en el cluster más grande.
library(corrplot)
par(mfrow = c(2,2))
for(nombre in names(resultados_sbm)){
corrplot::corrplot(resultados_sbm[[nombre]]$Pi, type = "full", col.lim = c(0,1), method = "shade",
addgrid.col = "gray90", tl.col = "black", title = nombre, mar = c(0,0,2,0))
}Loas visualizaciones de las probabilidades de interacción \(\pi_{q,r}\) muestran altos valores en las diagonales y valores cercanos a cero fuera de estas.
get_adjacency_ordered <- function(xi, A) {
xi2 <- xi[order(xi)]
indices <- order(xi)
d <- NULL
for (i in 1:(length(xi)-1)) if (xi2[i] != xi2[i+1]) d <- c(d, i)
list(A = A[indices,indices], d = d)
}
heat.plot0 <- function (mat, show.grid = FALSE, cex.axis, tick, labs, col.axis, ...){
JJ <- dim(mat)[1]
colorscale <- c("white", rev(heat.colors(100)))
if(missing(labs)) labs <- 1:JJ
if(missing(col.axis)) col.axis <- rep("black", JJ)
if(missing(cex.axis)) cex.axis <- 1
if(missing(tick)) tick <- TRUE
## adjacency matrix
image(seq(1, JJ), seq(1, JJ), mat, axes = FALSE, xlab = "", ylab = "", col = colorscale[seq(floor(100*min(mat)), floor(100*max(mat)))], ...)
for(j in 1:JJ){
axis(1, at = j, labels = labs[j], las = 2, cex.axis = cex.axis, tick, col.axis = col.axis[j], col.ticks = col.axis[j])
axis(2, at = j, labels = labs[j], las = 2, cex.axis = cex.axis, tick, col.axis = col.axis[j], col.ticks = col.axis[j])
}
box()
if(show.grid) grid(nx = JJ, ny = JJ)
}adj_ordenadas <- lapply(names(resultados_sbm), function(nombre){
get_adjacency_ordered(xi = resultados_sbm[[nombre]]$labs,
A = Y[[nombre]])
})
names(adj_ordenadas) <- names(resultados_sbm)
plot_sbm <- function(nombre){
tmp <- adj_ordenadas[[nombre]]
labs <- resultados_sbm[[nombre]]$labs
cols <- RColorBrewer::brewer.pal(5, "Set1")
par(mfrow = c(1,2), mar = c(1,1,1,1), oma = c(0,0,3,0))
set.seed(123)
plot(redes[[nombre]],
layout = layout_nicely,
vertex.label = NA,
vertex.size = 5,
vertex.color = cols[labs],
vertex.frame.color = cols[labs],
edge.color = adjustcolor("black",0.1))
heat.plot0(mat = tmp$A, tick = FALSE, labs = NA, asp = 1)
abline(v = tmp$d + .5, h = tmp$d + .5)
mtext(paste0("Agrupamiento de la Red de ", nombre, " - Modelo de bloques estocásticos"), side = 3, outer = TRUE, line = 1, cex = 1.2)
}
plot_sbm("uribe")Vemos que el modelo de bloques estocásticos agrupa a los funcionarios de acuerdo a su posición y conexiones en la red. Si la red de poder cuenta con varias díadas o grupos de díadas aisladas, estas las clasifica como un solo cluster, comportamiento bastante evidente en la red de Duque. Así mismo, si hay un grupo de funcionarios muy conectados entre sí (lo que se traduce a que han compartido varios cargos públicos entre ellos), el modelo los toma como un grupo. Del mismo modo, el modelo minimiza las conexiones entre grupos, como se observa en la matriz de adyacencia ordenada.
En todos los casos, este modelo da como resultado menos agrupaciones que los modelos de agrupamiento estudiados en el análisis estructural de cada red. Aunque, para el gobierno Santos esta distribución es bastante similar al algoritmo de comunidades.
Ahora, se calcula un modelo de grafos aleatorios exponenciales para cada una de las redes de poder. En estos, se incluyen loos efectos homofílicos de las variables nodales.
Para no generar ruido con la salida del modelo, se muestran solo los coeficientes estimados diferentes de -Inf.
library(ergm)
library(intergraph)
grafos <- lapply(redes, asNetwork)
set.seed(123)
ergm_modelos <- lapply(grafos, function(red){
ergm(red ~ edges + nodematch('PARTIDO_POLITICO', diff = T) + nodematch("GENERO") +
nodematch("CARRERA", diff = T) + nodematch("NIVEL_DE_ESTUDIOS", diff = T) ) # + absdiff('EDAD')
})
s <- summary(ergm_modelos$uribe)
as.data.frame(s$coefficients) %>%
filter(is.finite(Estimate)) %>%
mutate(Significativo = if_else(`Pr(>|z|)` < 0.05, "Sí", "No") )## Estimate Std. Error MCMC %
## edges -2.657171925 0.2127292 0
## nodematch.PARTIDO_POLITICO.CAMBIO RADICAL 17.390427266 559.1005420 0
## nodematch.PARTIDO_POLITICO.CENTRO DEMOCRATICO -0.008521488 0.4237105 0
## nodematch.PARTIDO_POLITICO.NA -0.271427163 0.2782563 0
## nodematch.PARTIDO_POLITICO.PARTIDO CONSERVADOR -0.292377855 0.6242832 0
## nodematch.GENERO 0.055047053 0.2129514 0
## nodematch.CARRERA.ABOGADO 3.324887981 1.2336218 0
## nodematch.CARRERA.DERECHO 1.253630655 0.3683960 0
## nodematch.CARRERA.ECONOMISTA 0.586319482 0.3087731 0
## nodematch.NIVEL_DE_ESTUDIOS.POSGRADO 0.021503364 0.2136586 0
## nodematch.NIVEL_DE_ESTUDIOS.PREGRADO -0.144818403 0.3945207 0
## z value Pr(>|z|)
## edges -12.49086321 8.374011e-36
## nodematch.PARTIDO_POLITICO.CAMBIO RADICAL 0.03110429 9.751864e-01
## nodematch.PARTIDO_POLITICO.CENTRO DEMOCRATICO -0.02011158 9.839544e-01
## nodematch.PARTIDO_POLITICO.NA -0.97545732 3.293335e-01
## nodematch.PARTIDO_POLITICO.PARTIDO CONSERVADOR -0.46834170 6.395403e-01
## nodematch.GENERO 0.25849581 7.960243e-01
## nodematch.CARRERA.ABOGADO 2.69522466 7.034119e-03
## nodematch.CARRERA.DERECHO 3.40294289 6.666421e-04
## nodematch.CARRERA.ECONOMISTA 1.89886810 5.758182e-02
## nodematch.NIVEL_DE_ESTUDIOS.POSGRADO 0.10064357 9.198334e-01
## nodematch.NIVEL_DE_ESTUDIOS.PREGRADO -0.36707428 7.135636e-01
## Significativo
## edges Sí
## nodematch.PARTIDO_POLITICO.CAMBIO RADICAL No
## nodematch.PARTIDO_POLITICO.CENTRO DEMOCRATICO No
## nodematch.PARTIDO_POLITICO.NA No
## nodematch.PARTIDO_POLITICO.PARTIDO CONSERVADOR No
## nodematch.GENERO No
## nodematch.CARRERA.ABOGADO Sí
## nodematch.CARRERA.DERECHO Sí
## nodematch.CARRERA.ECONOMISTA No
## nodematch.NIVEL_DE_ESTUDIOS.POSGRADO No
## nodematch.NIVEL_DE_ESTUDIOS.PREGRADO No
## [1] 805.8883
## attr(,"vcov")
## [1] 0
## [1] 864.53
## attr(,"vcov")
## [1] 0
s <- summary(ergm_modelos$santos)
as.data.frame(s$coefficients) %>%
filter(is.finite(Estimate)) %>%
mutate(Significativo = if_else(`Pr(>|z|)` < 0.05, "Sí", "No") )## Estimate Std. Error MCMC %
## edges -3.05903335 0.1726453 0
## nodematch.PARTIDO_POLITICO.CAMBIO RADICAL 3.09769465 0.6683213 0
## nodematch.PARTIDO_POLITICO.NA -0.25821253 0.2586877 0
## nodematch.PARTIDO_POLITICO.PARTIDO CONSERVADOR 0.34885090 0.2952189 0
## nodematch.PARTIDO_POLITICO.PARTIDO DE LA U 0.02023941 0.5459795 0
## nodematch.PARTIDO_POLITICO.PARTIDO LIBERAL 0.40006870 0.3491680 0
## nodematch.GENERO 0.55474193 0.1548257 0
## nodematch.CARRERA.ABOGADO 0.95209211 0.3674435 0
## nodematch.CARRERA.ADMINISTRADOR DE EMPRESAS 1.17132883 1.1329253 0
## nodematch.CARRERA.DERECHO 0.80103295 0.2220753 0
## nodematch.CARRERA.ECONOMISTA 0.73198443 0.1895028 0
## nodematch.CARRERA.INGENIERO CIVIL 1.83562371 1.2287123 0
## nodematch.CARRERA.INGENIERO INDUSTRIAL 0.71918720 1.0646491 0
## nodematch.NIVEL_DE_ESTUDIOS.POSGRADO -0.07405274 0.1422831 0
## nodematch.NIVEL_DE_ESTUDIOS.PREGRADO 0.55651464 0.2964465 0
## z value Pr(>|z|)
## edges -17.71860510 3.012910e-70
## nodematch.PARTIDO_POLITICO.CAMBIO RADICAL 4.63503833 3.568713e-06
## nodematch.PARTIDO_POLITICO.NA -0.99816315 3.182003e-01
## nodematch.PARTIDO_POLITICO.PARTIDO CONSERVADOR 1.18166872 2.373372e-01
## nodematch.PARTIDO_POLITICO.PARTIDO DE LA U 0.03706991 9.704293e-01
## nodematch.PARTIDO_POLITICO.PARTIDO LIBERAL 1.14577723 2.518873e-01
## nodematch.GENERO 3.58300931 3.396585e-04
## nodematch.CARRERA.ABOGADO 2.59112542 9.566262e-03
## nodematch.CARRERA.ADMINISTRADOR DE EMPRESAS 1.03389768 3.011840e-01
## nodematch.CARRERA.DERECHO 3.60703316 3.097181e-04
## nodematch.CARRERA.ECONOMISTA 3.86265795 1.121600e-04
## nodematch.CARRERA.INGENIERO CIVIL 1.49394101 1.351910e-01
## nodematch.CARRERA.INGENIERO INDUSTRIAL 0.67551571 4.993482e-01
## nodematch.NIVEL_DE_ESTUDIOS.POSGRADO -0.52046042 6.027427e-01
## nodematch.NIVEL_DE_ESTUDIOS.PREGRADO 1.87728550 6.047897e-02
## Significativo
## edges Sí
## nodematch.PARTIDO_POLITICO.CAMBIO RADICAL Sí
## nodematch.PARTIDO_POLITICO.NA No
## nodematch.PARTIDO_POLITICO.PARTIDO CONSERVADOR No
## nodematch.PARTIDO_POLITICO.PARTIDO DE LA U No
## nodematch.PARTIDO_POLITICO.PARTIDO LIBERAL No
## nodematch.GENERO Sí
## nodematch.CARRERA.ABOGADO Sí
## nodematch.CARRERA.ADMINISTRADOR DE EMPRESAS No
## nodematch.CARRERA.DERECHO Sí
## nodematch.CARRERA.ECONOMISTA Sí
## nodematch.CARRERA.INGENIERO CIVIL No
## nodematch.CARRERA.INGENIERO INDUSTRIAL No
## nodematch.NIVEL_DE_ESTUDIOS.POSGRADO No
## nodematch.NIVEL_DE_ESTUDIOS.PREGRADO No
## [1] 1909.403
## attr(,"vcov")
## [1] 0
## [1] 2002.078
## attr(,"vcov")
## [1] 0
s <- summary(ergm_modelos$duque)
as.data.frame(s$coefficients) %>%
filter(is.finite(Estimate)) %>%
mutate(Significativo = if_else(`Pr(>|z|)` < 0.05, "Sí", "No") )## Estimate Std. Error MCMC %
## edges -3.0830651 0.2873155 0
## nodematch.PARTIDO_POLITICO.CAMBIO RADICAL 0.7084691 1.1320766 0
## nodematch.PARTIDO_POLITICO.CENTRO DEMOCRATICO -0.0603425 0.4455709 0
## nodematch.PARTIDO_POLITICO.NA -0.1230504 0.3481779 0
## nodematch.GENERO -0.1514626 0.2408031 0
## nodematch.CARRERA.ABOGADO 2.2859598 0.6102038 0
## nodematch.CARRERA.COMUNICADOR SOCIAL 15.9134670 882.7435075 0
## nodematch.CARRERA.DERECHO 0.6579104 0.3780659 0
## nodematch.CARRERA.ECONOMISTA 0.8463700 0.4241927 0
## nodematch.CARRERA.MEDICO 18.0570071 882.7434461 0
## nodematch.NIVEL_DE_ESTUDIOS.POSGRADO -0.1333613 0.2978838 0
## nodematch.NIVEL_DE_ESTUDIOS.PREGRADO 1.8871285 0.5292999 0
## z value Pr(>|z|)
## edges -10.73058991 7.312784e-27
## nodematch.PARTIDO_POLITICO.CAMBIO RADICAL 0.62581375 5.314371e-01
## nodematch.PARTIDO_POLITICO.CENTRO DEMOCRATICO -0.13542740 8.922740e-01
## nodematch.PARTIDO_POLITICO.NA -0.35341231 7.237794e-01
## nodematch.GENERO -0.62898940 5.293560e-01
## nodematch.CARRERA.ABOGADO 3.74622321 1.795169e-04
## nodematch.CARRERA.COMUNICADOR SOCIAL 0.01802728 9.856171e-01
## nodematch.CARRERA.DERECHO 1.74020039 8.182384e-02
## nodematch.CARRERA.ECONOMISTA 1.99524912 4.601572e-02
## nodematch.CARRERA.MEDICO 0.02045556 9.836800e-01
## nodematch.NIVEL_DE_ESTUDIOS.POSGRADO -0.44769568 6.543728e-01
## nodematch.NIVEL_DE_ESTUDIOS.PREGRADO 3.56532970 3.633991e-04
## Significativo
## edges Sí
## nodematch.PARTIDO_POLITICO.CAMBIO RADICAL No
## nodematch.PARTIDO_POLITICO.CENTRO DEMOCRATICO No
## nodematch.PARTIDO_POLITICO.NA No
## nodematch.GENERO No
## nodematch.CARRERA.ABOGADO Sí
## nodematch.CARRERA.COMUNICADOR SOCIAL No
## nodematch.CARRERA.DERECHO No
## nodematch.CARRERA.ECONOMISTA Sí
## nodematch.CARRERA.MEDICO No
## nodematch.NIVEL_DE_ESTUDIOS.POSGRADO No
## nodematch.NIVEL_DE_ESTUDIOS.PREGRADO Sí
## [1] 604.6512
## attr(,"vcov")
## [1] 0
## [1] 669.259
## attr(,"vcov")
## [1] 0
s <- summary(ergm_modelos$petro)
as.data.frame(s$coefficients) %>%
filter(is.finite(Estimate)) %>%
mutate(Significativo = if_else(`Pr(>|z|)` < 0.05, "Sí", "No") )## Estimate Std. Error MCMC %
## edges -3.48533450 0.1669855 0
## nodematch.PARTIDO_POLITICO.COLOMBIA HUMANA -0.39970984 1.0555848 0
## nodematch.PARTIDO_POLITICO.NA 0.08513185 0.2316415 0
## nodematch.PARTIDO_POLITICO.PACTO HISTORICO -0.92643178 0.5916804 0
## nodematch.PARTIDO_POLITICO.PARTIDO DE LA U 2.03756212 1.2711793 0
## nodematch.PARTIDO_POLITICO.PARTIDO LIBERAL 1.54019457 0.5700816 0
## nodematch.GENERO 0.64856190 0.1670440 0
## nodematch.CARRERA.ABOGADO 1.17791139 0.3657250 0
## nodematch.CARRERA.COMUNICADOR SOCIAL 2.44829562 1.2604505 0
## nodematch.CARRERA.DERECHO 0.56832797 0.4118135 0
## nodematch.CARRERA.ECONOMISTA 0.49815883 0.4758046 0
## nodematch.CARRERA.INGENIERO CIVIL 2.42817242 1.2450279 0
## nodematch.CARRERA.MEDICO 16.87130033 324.7442605 0
## nodematch.CARRERA.POLITOLOGO 0.84151372 0.7627065 0
## nodematch.NIVEL_DE_ESTUDIOS.POSGRADO 0.10652895 0.1769714 0
## nodematch.NIVEL_DE_ESTUDIOS.PREGRADO 0.29463417 0.2402725 0
## z value Pr(>|z|)
## edges -20.87207354 9.607473e-97
## nodematch.PARTIDO_POLITICO.COLOMBIA HUMANA -0.37866198 7.049389e-01
## nodematch.PARTIDO_POLITICO.NA 0.36751547 7.132346e-01
## nodematch.PARTIDO_POLITICO.PACTO HISTORICO -1.56576385 1.174039e-01
## nodematch.PARTIDO_POLITICO.PARTIDO DE LA U 1.60289119 1.089587e-01
## nodematch.PARTIDO_POLITICO.PARTIDO LIBERAL 2.70170896 6.898412e-03
## nodematch.GENERO 3.88258098 1.033536e-04
## nodematch.CARRERA.ABOGADO 3.22075688 1.278526e-03
## nodematch.CARRERA.COMUNICADOR SOCIAL 1.94239735 5.208902e-02
## nodematch.CARRERA.DERECHO 1.38006154 1.675677e-01
## nodematch.CARRERA.ECONOMISTA 1.04698188 2.951079e-01
## nodematch.CARRERA.INGENIERO CIVIL 1.95029561 5.114090e-02
## nodematch.CARRERA.MEDICO 0.05195257 9.585665e-01
## nodematch.CARRERA.POLITOLOGO 1.10332578 2.698857e-01
## nodematch.NIVEL_DE_ESTUDIOS.POSGRADO 0.60195570 5.472036e-01
## nodematch.NIVEL_DE_ESTUDIOS.PREGRADO 1.22625002 2.201046e-01
## Significativo
## edges Sí
## nodematch.PARTIDO_POLITICO.COLOMBIA HUMANA No
## nodematch.PARTIDO_POLITICO.NA No
## nodematch.PARTIDO_POLITICO.PACTO HISTORICO No
## nodematch.PARTIDO_POLITICO.PARTIDO DE LA U No
## nodematch.PARTIDO_POLITICO.PARTIDO LIBERAL Sí
## nodematch.GENERO Sí
## nodematch.CARRERA.ABOGADO Sí
## nodematch.CARRERA.COMUNICADOR SOCIAL No
## nodematch.CARRERA.DERECHO No
## nodematch.CARRERA.ECONOMISTA No
## nodematch.CARRERA.INGENIERO CIVIL No
## nodematch.CARRERA.MEDICO No
## nodematch.CARRERA.POLITOLOGO No
## nodematch.NIVEL_DE_ESTUDIOS.POSGRADO No
## nodematch.NIVEL_DE_ESTUDIOS.PREGRADO No
## [1] 1339.126
## attr(,"vcov")
## [1] 0
## [1] 1437.108
## attr(,"vcov")
## [1] 0
Empezamos comparando las estimaciones de edges para las
redes de poder. Para la red de Uribe, la probabilidad base de que exista
una arista entre dos funcionarios que no comparten ninguna de las
características modeladas es de 6.55%, esta es de tan solo 4.48% para la
red de Santos, 4.48% para la red de Duque y 2.97% para la de Petro.
Pasando al efecto del partido político manteniendo lo demás constante observamos que, por ejemplo, dos funcionarios que hicieron parte de Cambio radical durante los gobiernos de Santos tuvieron una probabilidad estimada de conexión de 50.96% frente al 4.5% de un par que no compartía esta afiliación. Aunque se obtienen estimaciones de la homofilia de este partido también para los gobiernos de Uribe y Duque, estos no son significativos. En cuanto a la red del gobierno de Petro, dos funcionarios del Partido liberal tienen una probabilidad de conexión del 12.51%. Lo que nos dice que si dos personas hacen parte de este partido, la probabilidad de que ocupen los mismos cargos en este gobierno aumenta un 9.53%.
Continuamos con las carreras universitarias. Durante los gobiernos de Uribe, si dos funcionarios eran abogados la probabilidad de que ocuparan los mismos cargos era de 66.09% frente al 6.55% de aquellos que no comparten esta profesión. En esta misma red, estudiar derecho también tiene un efecto significativo: aumenta la probabilidad de conexión en un 4.64%. Algo similar ocurre en los gobiernos de Santos: la probabilidad de conexión aumenta (respecto de la probabilidad base) si dos funcionarios son abogados, economistas o estudiaron derecho. En general, siempre que dos funcionarios son abogados, es más probable que ocupen los mismos puestos públicos sin importar el gobierno de turno.
Las redes de Santos y Petro son las únicas en las que tener el mismo género aumenta la posibilidad de conexión de forma significativa.
Finalmente, en la red de Duque, tener un nivel de estudios de Pregrado aumenta la probabilidad de que dos funcionarios compartan cargos públicos. Esta es la única red que muestra un efecto homofilico significativo en el nivel de estudios.
## Analysis of Deviance Table
##
## Model 1: red ~ edges + nodematch("PARTIDO_POLITICO", diff = T) + nodematch("GENERO") +
## nodematch("CARRERA", diff = T) + nodematch("NIVEL_DE_ESTUDIOS",
## diff = T)
## Df Deviance Resid. Df Resid. Dev Pr(>|Chisq|)
## NULL 1527 2134.89
## Model 1: 11 1351 1516 783.89 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Analysis of Deviance Table
##
## Model 1: red ~ edges + nodematch("PARTIDO_POLITICO", diff = T) + nodematch("GENERO") +
## nodematch("CARRERA", diff = T) + nodematch("NIVEL_DE_ESTUDIOS",
## diff = T)
## Df Deviance Resid. Df Resid. Dev Pr(>|Chisq|)
## NULL 3563 4949.1
## Model 1: 15 3069.7 3548 1879.4 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Analysis of Deviance Table
##
## Model 1: red ~ edges + nodematch("PARTIDO_POLITICO", diff = T) + nodematch("GENERO") +
## nodematch("CARRERA", diff = T) + nodematch("NIVEL_DE_ESTUDIOS",
## diff = T)
## Df Deviance Resid. Df Resid. Dev Pr(>|Chisq|)
## NULL 1610 2291.54
## Model 1: 12 1710.9 1598 580.65 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Analysis of Deviance Table
##
## Model 1: red ~ edges + nodematch("PARTIDO_POLITICO", diff = T) + nodematch("GENERO") +
## nodematch("CARRERA", diff = T) + nodematch("NIVEL_DE_ESTUDIOS",
## diff = T)
## Df Deviance Resid. Df Resid. Dev Pr(>|Chisq|)
## NULL 3374 4717.6
## Model 1: 16 3410.4 3358 1307.1 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Para todas las redes de poder, haber incorporado la homofilia de las variables nodales implica una mejora en el ajuste respecto de un modelo nulo. Vemos que todos los modelos ERGM propuestos son significativos.
Por último, se ajusta un modelo de sociabilidad con variables nodales que representan efectos homofílicos. Estas variables se incorporaron en el modelo como \(x_{ijk}=I(\text{categoria}_i=\text{categoria}_j)\), donde \(k=1, ..., p\) es el ínidce de la variable nodal estudiada.
Primero, se muestran las distribuciones condicionales completas necesarias para el Muestreador de Gibbs.
# Distribuciones condicionales completas
# DCC 1: Muestreo de z
sample_z <- function(y, mu, delta, z, beta, X) {
n <- nrow(y)
for (i in 1:(n - 1)) {
for (j in (i + 1):n) {
mean_z <- mu + delta[i] + delta[j] + sum(X[i,j,]*beta)
if (y[i, j] == 1) {
z[i, j] <- truncnorm::rtruncnorm(n = 1, a = 0, b = Inf, mean = mean_z, sd = 1)
} else {
z[i, j] <- truncnorm::rtruncnorm(n = 1, a = -Inf, b = 0, mean = mean_z, sd = 1)
}
z[j, i] <- z[i, j] # Simetría
}
}
return(z)
}
# DCC 2: Muestreo de mu
sample_mu <- function(z, delta, sigma2, beta, X) {
x_beta <- apply(X, c(1,2), crossprod, beta)
indices_sup <- upper.tri(z)
v2_mu <- 1 / (1 / sigma2 + sum(indices_sup))
m_mu <- v2_mu * sum(z[indices_sup] - delta[row(z)[indices_sup]] - delta[col(z)[indices_sup]]- x_beta[indices_sup])
return(rnorm(1, mean = m_mu, sd = sqrt(v2_mu)))
}
# DCC 3: Muestreo de delta
sample_delta <- function(z, mu, tau2, delta, beta, X) {
n <- length(delta)
for (i in 1:n) {
neighbors <- setdiff(1:n, i)
v2_delta <- 1 / (1 / tau2 + length(neighbors))
x_beta_neighbors <- sapply(neighbors, function(j) sum(X[i, j, ] * beta))
m_delta <- v2_delta * sum(z[i, neighbors] - mu - delta[neighbors]- x_beta_neighbors)
delta[i] <- rnorm(1, mean = m_delta, sd = sqrt(v2_delta))
}
return(delta)
}
# DCC 5: Muestreo de beta
sample_beta <- function(z, mu, delta, X, mu0, Sigma0_inv) {
n <- nrow(z)
p <- dim(X)[3]
# Acumuladores para la FCD Normal Multivariada
sum_XXt <- matrix(0, p, p)
sum_Xw <- matrix(0, p, 1)
for (i in 1:(n - 1)) {
for (j in (i + 1):n) {
w_ij <- z[i, j] - mu - delta[i] - delta[j]
x_ij <- matrix(X[i,j,], ncol = 1)
sum_XXt <- sum_XXt + x_ij %*% t(x_ij)
sum_Xw <- sum_Xw + x_ij * w_ij
}
}
V_beta <- solve(Sigma0_inv + sum_XXt)
M_beta <- V_beta %*% (Sigma0_inv %*% mu0 + sum_Xw)
return(as.numeric(MASS::mvrnorm(1, mu = M_beta, Sigma = V_beta)))
}
# DCC 5: Muestreo de sigma^2
sample_sigma2 <- function(mu) {
a_sigma_post <- a_sigma + 0.5
b_sigma_post <- b_sigma + 0.5 * mu^2
return(1 / rgamma(1, shape = a_sigma_post, rate = b_sigma_post))
}
# DCC 6: Muestreo de tau^2
sample_tau2 <- function(delta) {
n <- length(delta)
a_tau_post <- a_tau + n / 2
b_tau_post <- b_tau + 0.5 * sum(delta^2)
return(1 / rgamma(1, shape = a_tau_post, rate = b_tau_post))
}# Muestreador de Gibbs
gibbs_sampler <- function(y, X, n_iter, n_burn, n_thin) {
n <- nrow(y)
p <- 4 # Número de covariables
# Inicialización
mu <- 0
delta <- rnorm(n, 0, 1)
beta <- rep(0, p)
sigma2 <- 1
tau2 <- 1
z <- matrix(0, n, n) # Variables auxiliares
# Almacenamiento
n_samples <- (n_iter - n_burn) / n_thin
samples <- list(mu = numeric(n_samples),
delta = matrix(0, nrow = n_samples, ncol = n),
beta = matrix(0, nrow = n_samples, ncol = p), #
sigma2 = numeric(n_samples),
tau2 = numeric(n_samples))
# Muestreo
cat("Iniciando muestreador de Gibbs...\n")
for (t in 1:n_iter) {
# Llamar a las funciones de muestreo
z <- sample_z(y, mu, delta, z, beta, X)
mu <- sample_mu(z, delta, sigma2, beta, X)
delta <- sample_delta(z, mu, tau2, delta, beta, X)
beta <- sample_beta(z, mu, delta, X, mu0, Sigma0_inv)
sigma2 <- sample_sigma2(mu)
tau2 <- sample_tau2(delta)
# Almacenar muestras según n_thin
if (t > n_burn && (t - n_burn) %% n_thin == 0) {
idx <- (t - n_burn) / n_thin
samples$mu[idx] <- mu
samples$delta[idx, ] <- delta
samples$beta[idx, ] <- beta
samples$sigma2[idx] <- sigma2
samples$tau2[idx] <- tau2
}
# Mostrar progreso
if (t %% (n_iter / 10) == 0) {
cat(sprintf("Progreso: %d%% completado\n", (t / n_iter) * 100))
}
}
cat("Muestreador completado.\n")
return(samples)
}df_atributos_total <- map_df(redes, function(g) {
# Extrae los nodos y sus atributos como un dataframe
igraph::as_data_frame(g, what = "vertices")
}, .id = "GOBIERNO")
atributos_uribe <- df_atributos_total %>%
filter(GOBIERNO == "uribe")
atributos_santos <- df_atributos_total %>%
filter(GOBIERNO == "santos")
atributos_duque <- df_atributos_total %>%
filter(GOBIERNO == "duque")
atributos_petro <- df_atributos_total %>%
filter(GOBIERNO == "petro")
# Función para obtener las variables homofílicas
homofilia <- function(atributo){
outer(atributo, atributo, FUN = "==") * 1
}
X_uribe <- array(0, dim = c(dim(Y_uribe)[1], dim(Y_uribe)[1], 4))
X_uribe[,,1] <- homofilia(atributos_uribe$GENERO)
X_uribe[,,2] <- homofilia(atributos_uribe$PARTIDO_POLITICO)
X_uribe[,,3] <- homofilia(atributos_uribe$NIVEL_DE_ESTUDIOS)
X_uribe[,,4] <- homofilia(atributos_uribe$CARRERA)
X_santos <- array(0, dim = c(dim(Y_santos)[1], dim(Y_santos)[1], 4))
X_santos[,,1] <- homofilia(atributos_santos$GENERO)
X_santos[,,2] <- homofilia(atributos_santos$PARTIDO_POLITICO)
X_santos[,,3] <- homofilia(atributos_santos$NIVEL_DE_ESTUDIOS)
X_santos[,,4] <- homofilia(atributos_santos$CARRERA)
X_duque <- array(0, dim = c(dim(Y_duque)[1], dim(Y_duque)[1], 4))
X_duque[,,1] <- homofilia(atributos_duque$GENERO)
X_duque[,,2] <- homofilia(atributos_duque$PARTIDO_POLITICO)
X_duque[,,3] <- homofilia(atributos_duque$NIVEL_DE_ESTUDIOS)
X_duque[,,4] <- homofilia(atributos_duque$CARRERA)
X_petro <- array(0, dim = c(dim(Y_petro)[1], dim(Y_petro)[1], 4))
X_petro[,,1] <- homofilia(atributos_petro$GENERO)
X_petro[,,2] <- homofilia(atributos_petro$PARTIDO_POLITICO)
X_petro[,,3] <- homofilia(atributos_petro$NIVEL_DE_ESTUDIOS)
X_petro[,,4] <- homofilia(atributos_petro$CARRERA)# Hiperparametros
p <- 4
a_sigma <- 2
b_sigma <- 1
a_tau <- 2
b_tau <- 1
mu0 <- rep(0, p)
Sigma0_inv <- diag(1, p)
# Ajustar el modelo usando Gibbs
n_iter <- 100000 + 10000
n_burn <- 10000
n_thin <- 5
set.seed(123)
samples_uribe <- gibbs_sampler(Y$uribe, X_uribe, n_iter, n_burn, n_thin)
samples_santos <- gibbs_sampler(Y$santos, X_santos, n_iter, n_burn, n_thin)
samples_duque <- gibbs_sampler(Y$duque, X_duque, n_iter, n_burn, n_thin)
samples_petro <- gibbs_sampler(Y$petro , X_petro, n_iter, n_burn, n_thin)
save(samples_uribe, file = "samples_socio_uribe.RData")
save(samples_santos, file = "samples_socio_santos.RData")
save(samples_duque, file = "samples_socio_duque.RData")
save(samples_petro, file = "samples_socio_petro.RData") # Cargar las muestras calculadas previamente
load("Datos/samples_socio_uribe.RData")
load("Datos/samples_socio_santos.RData")
load("Datos/samples_socio_duque.RData")
load("Datos/samples_socio_petro.RData")
muestras <- list(uribe = samples_uribe,
santos = samples_santos,
duque = samples_duque,
petro = samples_petro)
Xs <- list(uribe = X_uribe,
santos = X_santos,
duque = X_duque,
petro = X_petro)# Función para calcular la log-verosimilitud
log_likelihood <- function(y, X, samples) {
n <- nrow(y) # Número de nodos
log_lik_samples <- numeric(length(samples$mu)) # Almacenar la log-verosimilitud para cada muestra
for (s in seq_along(samples$mu)) {
mu <- samples$mu[s]
delta <- samples$delta[s, ]
beta <- samples$beta[s, ]
log_lik <- 0
for (i in 1:(n - 1)) {
for (j in (i + 1):n) {
eta_ij <- mu + delta[i] + delta[j] + sum(X[i,j,]*beta) # Predictor lineal
p_ij <- pnorm(eta_ij) # Probabilidad del modelo (probit)
# Sumar la contribución del par (i, j) a la log-verosimilitud
log_lik <- log_lik + y[i, j] * log(p_ij + 1e-10) + (1 - y[i, j]) * log(1 - p_ij + 1e-10)
}
}
log_lik_samples[s] <- log_lik
}
return(log_lik_samples)
}Siempre que se lleva a cabo un MCMC es importante revisar la convergencia. Por esto, se muestra a continuación la log-verosimilitud del modelo para cada red de poder.
# Calcular la log-verosimilitud para cada gobierno
#log_liks <- Map(log_likelihood, Y, Xs, muestras)
#save(log_liks, file = "log_lik_gobiernos.RData")
load("Datos/log_lik_gobiernos.RData")
par(mfrow = c(2, 2))
plot(log_liks$uribe, main = "Log-verosimilitud Gobiernos de Uribe", xlab = "Iteración", ylab = "Log-verosimilitud", pch = 20, col = "cyan3")
plot(log_liks$santos, main = "Log-verosimilitud Gobiernos de Santos", xlab = "Iteración", ylab = "Log-verosimilitud", pch = 20, col = "cyan3")
plot(log_liks$duque, main = "Log-verosimilitud Gobierno de Duque", xlab = "Iteración", ylab = "Log-verosimilitud", pch = 20, col = "cyan3")
plot(log_liks$petro, main = "Log-verosimilitud Gobierno de Petro", xlab = "Iteración", ylab = "Log-verosimilitud", pch = 20, col = "cyan3")Los gráficos de las log-verosimilitudes nos dejan ver que el modelo si convergió.
# Cargar librerías necesarias
suppressMessages(suppressWarnings(library(coda)))
# Función para calcular MCSE
mc_ee <- function(samples) {
n_eff <- effectiveSize(samples) # Tamaño efectivo
sd_est <- sd(as.vector(samples)) # Desviación estándar
mcse <- abs(sd_est / sqrt(n_eff)) # Error estándar de Monte Carlo (MCSE)
return(mcse)
}
lapply(muestras, function(x){
c(mu = mc_ee(x$mu),
sigma2 = mc_ee(x$sigma2),
tau2 = mc_ee(x$tau2) )
})## $uribe
## mu.var1 sigma2.var1 tau2.var1
## 0.0052010425 0.0177965131 0.0006749026
##
## $santos
## mu.var1 sigma2.var1 tau2.var1
## 0.0034124914 0.0200855450 0.0002159902
##
## $duque
## mu.var1 sigma2.var1 tau2.var1
## 0.0068161101 0.0194330862 0.0005770181
##
## $petro
## mu.var1 sigma2.var1 tau2.var1
## 0.0028607139 0.0190650696 0.0002757202
Del mismo modo, los errores estándar de Monte Carlo para los parámetros \(\mu\), \(\sigma^2\) y \(\tau^2\) toman valore pequeños, lo queb nos habla de una alta precisión en la estimación.
Ajustado el modelo y obtenidas las 20.000 muestras podemos hacer inferencia sobre los parámetros estimados
# Función para preparar los deltas
crear_delta_df <- function(samples){
delta_mean <- colMeans(samples$delta)
delta_ci95 <- apply(samples$delta, 2, quantile, probs = c(0.025, 0.975))
data.frame(nodo = 1:length(delta_mean),
media = delta_mean,
li = delta_ci95[1,],
ls = delta_ci95[2,]) %>%
mutate(significativo = (li > 0 | ls < 0)) %>%
arrange(media) %>%
mutate(posicion = dplyr::row_number() )
}
delta_dfs <- lapply(muestras, crear_delta_df)
graficar_delta <- function(delta_df, titulo){
ggplot(delta_df,
aes(x = posicion,
y = media,
color = significativo) ) +
geom_errorbar(aes(ymin = li,
ymax = ls),
width = 0) +
geom_point(size = 2.5) +
geom_hline(yintercept = 0, linetype = 2, color = "red") +
geom_text(aes(label = nodo), nudge_y = 0.05, size = 3) +
scale_color_manual(values = c("FALSE" = "grey70",
"TRUE" = "blue") ) +
labs(title = titulo,
x = NULL,
y = expression(delta)) +
theme_minimal() +
theme(legend.position = "none",
plot.title = element_text(hjust = 0.5, face = "bold") )
}
par(mfrow = c(2, 2))
graficar_delta(delta_dfs$uribe, "Inferencia sobre los parámetros delta \n Gobierno Uribe")Los parámetros \(\delta_i\) representan la sociabilidad individual de cada individuo, lo que es su propensión a formar conexiones en la red. Note que hay leves diferencias en los nodos que cuentan con \(\delta_i\) significativamente superiores a 0. La red de poder de Petro es la que cuenta con la mayor cantidad de funcionarios altamente sociables, mientras que la red de poder de Duque tiene la menor cantidad de nodos altamente sociables con solo 4. Se cuenta como “altamente sociables” a los nodos con \(\delta_i\) cuyos intervalos de credibilidad no comprenden el cero.
Por otro lado, en los gobiernos de Uribe y Santos hay presencia de funcionarios con una sociabilidad significativamente baja.
Siguiendo el análisis de la inferencia, se obtienen las probabilidades de interaciión \(\theta_{ij}=\Phi(\mu+\delta_i+\delta_j+\boldsymbol{\beta}^\top\boldsymbol{x}_{ij})\)
# Función para calcular la matriz de probabilidades theta_ij
compute_theta <- function(samples, X) {
n_samples <- length(samples$mu)
n <- ncol(samples$delta)
theta_avg <- matrix(0, n, n) # Inicializar matriz promedio
# Iterar sobre cada muestra
for (s in 1:n_samples) {
mu <- samples$mu[s]
delta <- samples$delta[s, ]
beta <- samples$beta[s, ]
theta <- matrix(0, n, n)
# Calcular theta_ij para cada par (i, j)
for (i in 1:(n - 1)) {
for (j in (i + 1):n) {
theta[i, j] <- pnorm(mu + delta[i] + delta[j] + sum(X[i,j,]*beta) ) # Probit link
theta[j, i] <- theta[i, j] # Simetría
}
}
# Promediar sobre las muestras
theta_avg <- theta_avg + theta/n_samples
}
return(theta_avg)
}
# Calcular la matriz promedio de theta_ij
#theta_avg <- Map(compute_theta, muestras, Xs)
#save(theta_avg, file = "theta_avg_gobiernos.RData")
load("Datos/theta_avg_gobiernos.RData")##
## Adjuntando el paquete: 'reshape2'
## The following object is masked from 'package:tidyr':
##
## smiths
library(patchwork)
graficar_theta_vs_y <- function(theta, Y, titulo = "") {
theta_df <- reshape2::melt(theta)
names(theta_df) <- c("Nodo_i", "Nodo_j", "Valor")
y_df <- reshape2::melt(Y)
names(y_df) <- c("Nodo_i", "Nodo_j", "Valor")
# Heatmap probabilidades posteriores
g_theta <- ggplot(theta_df,
aes(x = Nodo_i,
y = Nodo_j,
fill = Valor)) +
geom_tile() +
coord_equal() +
scale_fill_gradient(low = "white",
high = "blueviolet") +
labs(title = paste("Probabilidades de interacción - ", titulo),
x = "Nodo i",
y = "Nodo j") +
theme_minimal()
# Heatmap matriz observada
g_y <- ggplot(y_df,
aes(x = Nodo_i,
y = Nodo_j,
fill = factor(Valor)) ) +
geom_tile() +
coord_equal() +
scale_fill_manual(values = c("0" = "white",
"1" = "blue") ) +
labs(title = paste("Matriz observada - ", titulo),
x = "Nodo i",
y = "Nodo j") +
theme_minimal() +
theme(
axis.text.x = element_blank(),
axis.text.y = element_blank(),
axis.ticks = element_blank(),
legend.position = "none" )
g_theta + g_y
}
graficos_theta <- Map(graficar_theta_vs_y, theta_avg, Y, names(theta_avg) )
graficos_theta$uribeLas matrices de probabilidades de interacción construidas por el modelo de sociabilidad muestran un ajuste a las matrices de adyacencia observadas. Las probabilidades de interacción parecen tomar valores más altos para los \((ij)\) que dieron cuenta de la existencia de un enlace, sin asignar una probabilidad de cero a las demás zonas; mostrando que no habría un sobre ajuste del modelo a los datos.
Una vez ajustados los modelos, es momento de evaluar la bondad de ajuste de estos.
# - - - - - - Función para Bondad de ajuste (SBM) - - - - - - - -
gof_sbm <- function(red, modelo_sbm, n_sim = 1000){
# Número óptimo de bloques
Q <- which.max(modelo_sbm$ICL)
# Pertenencias
Z <- modelo_sbm$memberships[[Q]]$Z
labs <- apply(Z, 1, which.max)
# Parámetros SBM
alpha <- table(labs) / vcount(red)
Pi <- modelo_sbm$model_parameters[[Q]]$pi
Pi <- (Pi + t(Pi))/2
# Estadísticas observadas
observed <- c(Asort_Genero = assortativity_nominal(red,
types = as.integer(factor(vertex_attr(red, "GENERO"))),
directed = FALSE),
Asort_Partido = assortativity_nominal(red,
types = as.integer(factor(vertex_attr(red, "PARTIDO_POLITICO"))),
directed = FALSE),
Asort_Estudios = assortativity_nominal(red,
types = as.integer(factor(vertex_attr(red, "NIVEL_DE_ESTUDIOS"))),
directed = FALSE),
Asort_Carrera = assortativity_nominal(red,
types = as.integer(factor(vertex_attr(red, "CARRERA"))),
directed = FALSE),
Densidad = edge_density(red),
Dist_Geodesica = mean_distance(red),
Transitividad = transitivity(red, type = "global"))
# Simulaciones
sim_stats <- data.frame(Asort_Genero = numeric(n_sim),
Asort_Partido = numeric(n_sim),
Asort_Estudios = numeric(n_sim),
Asort_Carrera = numeric(n_sim),
Densidad = numeric(n_sim),
Dist_Geodesica = numeric(n_sim),
Transitividad = numeric(n_sim) )
colnames(sim_stats) <- names(observed)
set.seed(123)
for(i in seq_len(n_sim)){
bs <- rmultinom(n = 1, size = vcount(red), prob = alpha)
g_sim <- igraph::sample_sbm(n = vcount(red), pref.matrix = Pi, block.sizes = bs, directed = FALSE)
# Estadísticas simuladas
sim_stats$Asort_Genero[i] <- assortativity_nominal(g_sim,
types = as.integer(factor(vertex_attr(red, "GENERO"))),
directed = FALSE)
sim_stats$Asort_Partido[i] <- assortativity_nominal(g_sim,
types = as.integer(factor(vertex_attr(red, "PARTIDO_POLITICO"))),
directed = FALSE)
sim_stats$Asort_Estudios[i] <- assortativity_nominal(g_sim,
types = as.integer(factor(vertex_attr(red, "NIVEL_DE_ESTUDIOS"))),
directed = FALSE)
sim_stats$Asort_Carrera[i] <- assortativity_nominal(g_sim,
types = as.integer(factor(vertex_attr(red, "CARRERA"))),
directed = FALSE)
sim_stats$Densidad[i] <- edge_density(g_sim)
sim_stats$Dist_Geodesica[i] <- mean_distance(g_sim)
sim_stats$Transitividad[i] <- transitivity(g_sim, type = "global")
}
#sim_stats <- as.data.frame(sim_stats)
# Tabla resumen
tabla <- data.frame(Estadistica = names(observed),
Observado = observed,
Media_Simulada = sapply(sim_stats, mean, na.rm = TRUE),
SD_Simulada = sapply(sim_stats, sd, na.rm = TRUE),
IC_Inf = sapply(sim_stats, quantile, probs = .025, na.rm = TRUE),
IC_Sup = sapply(sim_stats, quantile, probs = .975, na.rm = TRUE))
# Datos para gráfico
plot_df <- sim_stats %>%
pivot_longer(everything(), names_to = "Estadistica", values_to = "Valor")
obs_df <- data.frame(Estadistica = names(observed),
Observado = observed)
ci_df <- plot_df %>%
group_by(Estadistica) %>%
summarise(IC_Inf = quantile(Valor, .025),
IC_Sup = quantile(Valor, .975))
grafico <- ggplot(plot_df,
aes(x = Valor)) +
geom_histogram(#aes(y = after_stat(density)),
bins = 20,
fill = "grey85", color = "grey60") +
labs(x = "Valor", y = "Frecuencia") +
geom_vline(data = obs_df,
aes(xintercept = Observado),
color = "red", linewidth = 1) +
geom_vline(data = ci_df,
aes(xintercept = IC_Inf),
colour = "blue",
linetype = "dashed", linewidth = .8) +
geom_vline(data = ci_df,aes(xintercept = IC_Sup),
colour = "blue",
linetype = "dashed", linewidth = .8) +
facet_wrap(~Estadistica, scales = "free", ncol = 3) +
theme_minimal()
tabla_kable <- tabla %>%
kbl(caption = "Resumen de las estadísticas observadas y simuladas con IC del 95%", digits = 4,
row.names = FALSE,
col.names = c("Estadística", "Observado", "Media", "Desv. Est.", "IC Inferior", "IC Superior")
) %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE) #%>%
#column_spec(1, bold = TRUE)
list(tabla = tabla_kable,
grafico = grafico,
simulaciones = sim_stats)
}
gof_modelos_sbm <- mapply(FUN = gof_sbm,
red = redes,
modelo_sbm = sbm_modelos,
SIMPLIFY = FALSE)for(gob in names(gof_modelos_sbm)){
print(gof_modelos_sbm[[gob]]$grafico +
ggtitle(paste("Bondad de ajuste SBM -", gob)))
}for(gob in names(gof_modelos_sbm)){
cat("\n\n")
cat("Análisis bondad de ajuste - ", toupper(gob), "\n")
print(gof_modelos_sbm[[gob]]$tabla)
}| Estadística | Observado | Media | Desv. Est. | IC Inferior | IC Superior |
|---|---|---|---|---|---|
| Asort_Genero | 0.1396 | 0.0179 | 0.0757 | -0.1364 | 0.1723 |
| Asort_Partido | 0.0593 | -0.0295 | 0.0433 | -0.1131 | 0.0597 |
| Asort_Estudios | -0.0689 | -0.0532 | 0.0638 | -0.1793 | 0.0702 |
| Asort_Carrera | 0.1025 | -0.0024 | 0.0299 | -0.0624 | 0.0564 |
| Densidad | 0.0753 | 0.0785 | 0.0148 | 0.0545 | 0.1117 |
| Dist_Geodesica | 2.8144 | 3.6774 | 0.5745 | 2.7643 | 5.0720 |
| Transitividad | 0.6295 | 0.4064 | 0.0811 | 0.2430 | 0.5573 |
| Estadística | Observado | Media | Desv. Est. | IC Inferior | IC Superior |
|---|---|---|---|---|---|
| Asort_Genero | 0.1361 | 0.0367 | 0.0579 | -0.0675 | 0.1572 |
| Asort_Partido | 0.0395 | -0.0035 | 0.0224 | -0.0465 | 0.0409 |
| Asort_Estudios | 0.0756 | 0.2129 | 0.0711 | 0.0615 | 0.3439 |
| Asort_Carrera | 0.1246 | 0.0690 | 0.0345 | 0.0021 | 0.1353 |
| Densidad | 0.0779 | 0.0808 | 0.0086 | 0.0669 | 0.1006 |
| Dist_Geodesica | 2.8729 | 3.0859 | 0.3205 | 2.5731 | 3.7375 |
| Transitividad | 0.6386 | 0.3833 | 0.0672 | 0.2665 | 0.5293 |
| Estadística | Observado | Media | Desv. Est. | IC Inferior | IC Superior |
|---|---|---|---|---|---|
| Asort_Genero | -0.0157 | 0.0175 | 0.1112 | -0.1943 | 0.2468 |
| Asort_Partido | 0.0096 | -0.0367 | 0.0568 | -0.1519 | 0.0735 |
| Asort_Estudios | 0.3134 | -0.0411 | 0.0991 | -0.1961 | 0.1796 |
| Asort_Carrera | 0.1373 | -0.0267 | 0.0468 | -0.1181 | 0.0664 |
| Densidad | 0.0466 | 0.0464 | 0.0082 | 0.0327 | 0.0635 |
| Dist_Geodesica | 2.8486 | 4.3615 | 0.8740 | 2.6832 | 6.2944 |
| Transitividad | 0.6290 | 0.1768 | 0.0602 | 0.0600 | 0.2868 |
| Estadística | Observado | Media | Desv. Est. | IC Inferior | IC Superior |
|---|---|---|---|---|---|
| Asort_Genero | 0.2152 | 0.0320 | 0.0717 | -0.1169 | 0.1658 |
| Asort_Partido | 0.0506 | -0.0229 | 0.0299 | -0.0826 | 0.0352 |
| Asort_Estudios | 0.0983 | 0.0599 | 0.0675 | -0.0716 | 0.1869 |
| Asort_Carrera | 0.0471 | 0.0156 | 0.0233 | -0.0315 | 0.0627 |
| Densidad | 0.0505 | 0.0517 | 0.0089 | 0.0382 | 0.0714 |
| Dist_Geodesica | 2.8897 | 3.7846 | 0.3760 | 3.1192 | 4.6117 |
| Transitividad | 0.6515 | 0.3983 | 0.1537 | 0.1304 | 0.6926 |
A nivel general, el modelo de bloques estocásticos captura la densidad correctamente para todas las redes de poder de los presidentes estudiados, ya que vemos que el valor observado (línea roja) se encuentra dentro de los límites de credibilidad al \(95\%\) de las densidades de las redes simuladas, además de ubicarse casi en el centro de dichas distribuciones. Las demás características se analizan a continuación para cada gobierno.
Red de poder de Uribe: Para la red de poder de Uribe este modelo no captura correctamente la asortatividad de la variable ‘Carrera Universitaria’ ni la transitividad, las cuales subestima. También tiende a subestimar la asortatividad de ‘Género’ y ‘Partido político’, aunque medio las caprura (están al borde de los límites de credibilidad). Lo que si logra capturar es la homofilia en el ‘Nivel de estudios’
Red de poder de Santos: Aparte de la densidad, para esta red el modelo solo captura la distancia geodésica promedio. En este caso, el modelo subestima la transitividad y la homofilia de las variables ‘Carrera’, ‘Genero’ y ‘Partido politico’.
Red de poder de Duque: El modelo subestima la transitividad y la asortatividad de las variables ‘Carrera’ y ‘Nivel de estudios’. A diferencia de las redes anteriores, en esta si se reproduce correctamente la asortatividad de las varisbles ‘Genero’ y ‘Partido politico’. Finalmente, el modelo sobre-estima la distancia geodésica promedio.
Red de poder de Petro: En esta ocasión, el modelo alcanza a capturar la transitividad aunque la sigue subestimando; subestima la asortatividad de las variables ‘Género’ y ‘Partido político’. La asortatividad de ‘Nivel de estudios’ si se ve bien reproducida.
Retomando los análisis de todas las redes de poder, el modelo de bloques estocásticos no reproduce bien características estructurales de los grafos como la transitividad y la distancia geodésica promedio. Así mismo, como era de esperarse, se queda corto en caracterizar los efectos homofílicos de las variables nodales ya que estas no se incluyeron en la estimación.
Complementariamente, note que la distancia geodésica promedio fue la característica estructural con mayor variabilidad en las redes simuladas.
gof_ergm <- function(red, modelo_ergm, n_sim = 1000){
# Estadísticas observadas
obs_stats <- c(Asort_Carrera = assortativity_nominal(red,
types = as.integer(factor(vertex_attr(red, "CARRERA"))),
directed = FALSE),
Asort_Estudios = assortativity_nominal(red,
types = as.integer(factor(vertex_attr(red, "NIVEL_DE_ESTUDIOS"))),
directed = FALSE),
Asort_Genero = assortativity_nominal(red,
types = as.integer(factor(vertex_attr(red, "GENERO"))),
directed = FALSE),
Asort_Partido = assortativity_nominal(red,
types = as.integer(factor(vertex_attr(red, "PARTIDO_POLITICO"))),
directed = FALSE),
Densidad = edge_density(red),
Dist_geo_prom = mean_distance(red, directed = FALSE),
Transitividad = transitivity(red, type = "global")
)
# Simular redes
sim_networks <- simulate(modelo_ergm, nsim = n_sim, output = "network")
sim_stats <- matrix(NA, nrow = n_sim, ncol = length(obs_stats))
colnames(sim_stats) <- names(obs_stats)
for(i in seq_len(n_sim)){
g_sim <- intergraph::asIgraph(sim_networks[[i]])
sim_stats[i, ] <- c(assortativity_nominal(g_sim,
types = as.integer(factor(vertex_attr(g_sim, "CARRERA"))),
directed = FALSE),
assortativity_nominal(g_sim,
types = as.integer(factor(vertex_attr(g_sim, "NIVEL_DE_ESTUDIOS"))),
directed = FALSE),
assortativity_nominal(g_sim,
types = as.integer(factor(vertex_attr(g_sim, "GENERO"))),
directed = FALSE),
assortativity_nominal(g_sim,
types = as.integer(factor(vertex_attr(g_sim, "PARTIDO_POLITICO"))),
directed = FALSE),
edge_density(g_sim),
mean_distance(g_sim, directed = FALSE),
transitivity(g_sim, type = "global")
)
}
# Resumen
sim_mean <- apply(sim_stats, 2, mean, na.rm = TRUE)
sim_sd <- apply(sim_stats, 2, sd, na.rm = TRUE)
sim_ci_low <- apply(sim_stats, 2, quantile, probs = .025, na.rm = TRUE)
sim_ci_high <- apply(sim_stats, 2, quantile, probs = .975, na.rm = TRUE)
summary_df <- data.frame(Observed = obs_stats,
SimMean = sim_mean,
SimSD = sim_sd,
CI_low = sim_ci_low,
CI_high = sim_ci_high)
# Tabla
tabla <- data.frame(Estadistica = names(obs_stats),
Observado = obs_stats,
'Media' = sim_mean,
'Desv. Est.' = sim_sd,
'IC Inferior' = sim_ci_low,
'IC Superior' = sim_ci_high)
# Datos para ggplot
plot_df <- as.data.frame(sim_stats) %>%
pivot_longer(everything(),
names_to = "Estadistica",
values_to = "Valor")
obs_df <- data.frame(Estadistica = names(obs_stats),
Observado = obs_stats)
ci_df <- data.frame(Estadistica = names(obs_stats),
CI_low = sim_ci_low,
CI_high = sim_ci_high)
grafico <- ggplot(plot_df,
aes(x = Valor)) +
geom_histogram(bins = 20,
#aes(y = after_stat(count / sum(count))),
fill = "grey85",
colour = "grey50") +
labs(x = "Valor", y = "Frecuencia") +
geom_vline(data = obs_df,
aes(xintercept = Observado),
colour = "red",
linewidth = 1) +
geom_vline(data = ci_df,
aes(xintercept = CI_low),
colour = "blue",
linewidth = 0.8,
linetype = "dashed") +
geom_vline(data = ci_df,
aes(xintercept = CI_high),
colour = "blue",
linewidth = 0.8,
linetype = "dashed") +
facet_wrap(~Estadistica, scales = "free", ncol = 3) +
theme_minimal()
list(tabla = tabla,
summary = summary_df,
simulations = sim_stats,
grafico = grafico)
}
gof_modelos_ergm <- mapply(gof_ergm, redes, ergm_modelos, SIMPLIFY = FALSE)
gof_modelos_ergm$uribe$grafico +
ggtitle(paste("Resumen de las estadísticas observadas y simuladas con IC del 95%"))gof_modelos_ergm$santos$grafico +
ggtitle(paste("Resumen de las estadísticas observadas y simuladas con IC del 95%"))gof_modelos_ergm$duque$grafico +
ggtitle(paste("Resumen de las estadísticas observadas y simuladas con IC del 95%"))gof_modelos_ergm$petro$grafico +
ggtitle(paste("Resumen de las estadísticas observadas y simuladas con IC del 95%"))for(gob in names(gof_modelos_ergm)){
cat("\n\n")
cat("Análisis bondad de ajuste - ", toupper(gob), "\n")
print(gof_modelos_ergm[[gob]]$tabla %>%
kbl(row.names = FALSE, digits = 4, caption = paste("Bondad de ajuste ERGM -", gob),
col.names = c("Estadística", "Observado", "Media", "Desv. Est.", "IC Inferior", "IC Superior")
) %>%
kable_styling(
bootstrap_options = c("striped","hover","condensed"), full_width = FALSE) )
}| Estadística | Observado | Media | Desv. Est. | IC Inferior | IC Superior |
|---|---|---|---|---|---|
| Asort_Carrera | 0.1025 | 0.0797 | 0.0379 | 0.0083 | 0.1564 |
| Asort_Estudios | -0.0689 | -0.0500 | 0.0719 | -0.1896 | 0.0918 |
| Asort_Genero | 0.1396 | -0.0294 | 0.0869 | -0.1905 | 0.1480 |
| Asort_Partido | 0.0593 | -0.0117 | 0.0456 | -0.0923 | 0.0778 |
| Densidad | 0.0753 | 0.0756 | 0.0063 | 0.0623 | 0.0877 |
| Dist_geo_prom | 2.8144 | 2.9047 | 0.1620 | 2.6591 | 3.2754 |
| Transitividad | 0.6295 | 0.0869 | 0.0213 | 0.0461 | 0.1291 |
| Estadística | Observado | Media | Desv. Est. | IC Inferior | IC Superior |
|---|---|---|---|---|---|
| Asort_Carrera | 0.1246 | 0.1175 | 0.0282 | 0.0608 | 0.1723 |
| Asort_Estudios | 0.0756 | 0.0749 | 0.0569 | -0.0409 | 0.1850 |
| Asort_Genero | 0.1361 | 0.1393 | 0.0678 | 0.0068 | 0.2714 |
| Asort_Partido | 0.0395 | 0.0251 | 0.0241 | -0.0216 | 0.0712 |
| Densidad | 0.0779 | 0.0780 | 0.0045 | 0.0692 | 0.0874 |
| Dist_geo_prom | 2.8729 | 2.5579 | 0.0672 | 2.4339 | 2.6972 |
| Transitividad | 0.6386 | 0.0868 | 0.0124 | 0.0624 | 0.1108 |
| Estadística | Observado | Media | Desv. Est. | IC Inferior | IC Superior |
|---|---|---|---|---|---|
| Asort_Carrera | 0.1373 | 0.1323 | 0.0523 | 0.0298 | 0.2401 |
| Asort_Estudios | 0.3134 | 0.3133 | 0.1205 | 0.0689 | 0.5392 |
| Asort_Genero | -0.0157 | -0.0774 | 0.1075 | -0.2763 | 0.1412 |
| Asort_Partido | 0.0096 | -0.0387 | 0.0575 | -0.1498 | 0.0717 |
| Densidad | 0.0466 | 0.0465 | 0.0051 | 0.0369 | 0.0569 |
| Dist_geo_prom | 2.8486 | 3.8605 | 0.3910 | 3.2753 | 4.7863 |
| Transitividad | 0.6290 | 0.0580 | 0.0301 | 0.0000 | 0.1208 |
| Estadística | Observado | Media | Desv. Est. | IC Inferior | IC Superior |
|---|---|---|---|---|---|
| Asort_Carrera | 0.0471 | 0.0600 | 0.0255 | 0.0146 | 0.1096 |
| Asort_Estudios | 0.0983 | 0.0704 | 0.0657 | -0.0616 | 0.1941 |
| Asort_Genero | 0.2152 | 0.2572 | 0.0786 | 0.0991 | 0.4051 |
| Asort_Partido | 0.0506 | -0.0060 | 0.0299 | -0.0642 | 0.0547 |
| Densidad | 0.0505 | 0.0505 | 0.0038 | 0.0432 | 0.0582 |
| Dist_geo_prom | 2.8897 | 3.2047 | 0.1499 | 2.9394 | 3.5380 |
| Transitividad | 0.6515 | 0.0536 | 0.0152 | 0.0234 | 0.0838 |
El modelo de grafos aleatorias exponencial muestra una mejora en la reproducción de la homofilia en la red, respecto del modelo de bloques estocásticos. Esto se debe a que en su formulación, estaban incluidos los efectos homofílicos de las variables nodales de la red.
La red de poder de Uribe es la única en la que la asortatividad de la variable ‘Genero’ se encuentra al borde del intervalo de credibilidad de las estaísticas simuladas, en las demás redes el modelo captura correctamente esta característica.
La distancia geodésica es una estadística que se comporta diferente para cada red. En el caso de Uribe, el modelo si la captura bien; mientras que para los gobiernos de Duque y Petro el modelo sobre-estima esta característica de distancia. COntrariamente a lo que se observa en la red de Santos: una distancia geodésica subestimada por el modelo.
Aunque, como ya se mencionó, se ven mejoras en el ajuste de este modelo, este está muy lejos de capturar lo altamente transitividas que son las redes de poder.
gof_sociabilidad <- function(samples, g_obs, X, n_sim = 100) {
# estadísticos observados
obs_stats <- c(Asort_Carrera = assortativity_nominal(g_obs,
types = as.integer(factor(vertex_attr(g_obs, "CARRERA"))),
directed = FALSE),
Asort_Estudios = assortativity_nominal(g_obs,
types = as.integer(factor(vertex_attr(g_obs, "NIVEL_DE_ESTUDIOS"))),
directed = FALSE),
Asort_Genero = assortativity_nominal(g_obs,
types = as.integer(factor(vertex_attr(g_obs, "GENERO"))),
directed = FALSE),
Asort_Partido = assortativity_nominal(g_obs,
types = as.integer(factor(vertex_attr(g_obs, "PARTIDO_POLITICO"))),
directed = FALSE),
Densidad = edge_density(g_obs),
Dist_geodesica = mean_distance(g_obs, directed = FALSE),
Transitividad = transitivity(g_obs, type = "global")
)
n <- vcount(g_obs)
sim_stats <- matrix(NA, nrow = n_sim, ncol = length(obs_stats))
colnames(sim_stats) <- names(obs_stats)
N <- length(samples$mu)
for (s in 1:n_sim) {
iter <- sample(1:N, 1)
mu <- samples$mu[iter]
delta <- samples$delta[iter, ]
beta <- samples$beta[iter, ]
P <- matrix(0, n, n)
for (i in 1:(n-1)) {
for (j in (i+1):n) {
P[i,j] <- pnorm(mu + delta[i] + delta[j] + sum(X[i,j,] * beta))
P[j,i] <- P[i,j]
}
}
diag(P) <- 0
Y_sim <- matrix(rbinom(n*n, 1, as.vector(P)), n, n)
Y_sim[lower.tri(Y_sim)] <- t(Y_sim)[lower.tri(Y_sim)] # pisa lower con upper
diag(Y_sim) <- 0
g_sim <- graph_from_adjacency_matrix(Y_sim, mode = "undirected")
sim_stats[s, ] <- c(assortativity_nominal(g_sim,
types = as.integer(factor(vertex_attr(g_obs, "CARRERA"))),
directed = FALSE),
assortativity_nominal(g_sim,
types = as.integer(factor(vertex_attr(g_obs, "NIVEL_DE_ESTUDIOS"))),
directed = FALSE),
assortativity_nominal(g_sim,
types = as.integer(factor(vertex_attr(g_obs, "GENERO"))),
directed = FALSE),
assortativity_nominal(g_sim,
types = as.integer(factor(vertex_attr(g_obs, "PARTIDO_POLITICO"))),
directed = FALSE),
edge_density(g_sim),
mean_distance(g_sim, directed = FALSE),
transitivity(g_sim, type = "global") )
}
# resumen
sim_mean <- apply(sim_stats, 2, mean, na.rm = TRUE)
sim_sd <- apply(sim_stats, 2, sd, na.rm = TRUE)
sim_ci_low <- apply(sim_stats, 2, quantile, 0.025, na.rm = TRUE)
sim_ci_high <- apply(sim_stats, 2, quantile, 0.975, na.rm = TRUE)
tabla <- data.frame(Estadistica = names(obs_stats),
Observado = obs_stats,
Media_Simulada = sim_mean,
SimSD = sim_sd,
IC_Inf = sim_ci_low,
IC_Sup = sim_ci_high )
tabla_kable <- tabla %>%
kbl(caption = "Resumen de las estadísticas observadas y simuladas con IC del 95%",
digits = 4, row.names = FALSE,
col.names = c("Estadística", "Observado", "Media", "Desv. Est.", "IC Inferior", "IC Superior") ) %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE) #%>%
#column_spec(1, bold = TRUE)
return(list(tabla = tabla,
tabla_kable = tabla_kable,
simulations = sim_stats)
)
}
gof_all <- Map(
function(samples, g_obs, X) {gof_sociabilidad(samples, g_obs, X, n_sim = 1000)},
muestras, redes, Xs)
plot_gof_gg <- function(gof_obj, name = "") {
df_sum <- gof_obj$tabla
sim <- as.data.frame(gof_obj$simulations)
sim_long <- sim %>%
mutate(sim_id = 1:n()) %>%
pivot_longer(-sim_id,
names_to = "Estadistica",
values_to = "Valor")
obs_df <- data.frame(
Estadistica = rownames(df_sum),
Observado = df_sum$Observado,
CI_low = df_sum$IC_Inf,
CI_high = df_sum$IC_Inf)
ggplot(sim_long, aes(x = Valor)) +
geom_histogram(
bins = 20,
#aes(y = after_stat(count / sum(count))),
fill = "grey85",
colour = "grey50") +
geom_vline(data = obs_df,
aes(xintercept = Observado),
colour = "red",
linewidth = 1) +
geom_vline(data = obs_df,
aes(xintercept = CI_low),
colour = "blue",
linetype = "dashed") +
geom_vline(data = obs_df,
aes(xintercept = CI_high),
colour = "blue",
linetype = "dashed") +
facet_wrap(~Estadistica, scales = "free", ncol = 3) +
labs(x = "Valor", y = "Frecuencia") +
theme_minimal()
}
plot_gof_gg(gof_all$uribe, name = "Bondad de ajuste sociabilidad - Red Uribe")for(gob in names(gof_all)){
cat("Análisis bondad de ajuste - ", toupper(gob), "\n")
print(gof_all[[gob]]$tabla_kable)
}| Estadística | Observado | Media | Desv. Est. | IC Inferior | IC Superior |
|---|---|---|---|---|---|
| Asort_Carrera | 0.1025 | 0.0844 | 0.0476 | -0.0037 | 0.1823 |
| Asort_Estudios | -0.0689 | -0.0721 | 0.0938 | -0.2531 | 0.1157 |
| Asort_Genero | 0.1396 | 0.0950 | 0.1191 | -0.1360 | 0.3258 |
| Asort_Partido | 0.0593 | 0.0497 | 0.0639 | -0.0752 | 0.1764 |
| Densidad | 0.0753 | 0.0770 | 0.0090 | 0.0597 | 0.0955 |
| Dist_geodesica | 2.8144 | 2.6754 | 0.1557 | 2.4230 | 3.0287 |
| Transitividad | 0.6295 | 0.1443 | 0.0325 | 0.0830 | 0.2079 |
| Estadística | Observado | Media | Desv. Est. | IC Inferior | IC Superior |
|---|---|---|---|---|---|
| Asort_Carrera | 0.1246 | 0.1182 | 0.0381 | 0.0441 | 0.1939 |
| Asort_Estudios | 0.0756 | 0.0362 | 0.0703 | -0.0962 | 0.1704 |
| Asort_Genero | 0.1361 | 0.1232 | 0.0867 | -0.0367 | 0.3095 |
| Asort_Partido | 0.0395 | 0.0263 | 0.0343 | -0.0372 | 0.0960 |
| Densidad | 0.0779 | 0.0792 | 0.0060 | 0.0675 | 0.0916 |
| Dist_geodesica | 2.8729 | 2.5127 | 0.0810 | 2.3732 | 2.6847 |
| Transitividad | 0.6386 | 0.1193 | 0.0166 | 0.0876 | 0.1536 |
| Estadística | Observado | Media | Desv. Est. | IC Inferior | IC Superior |
|---|---|---|---|---|---|
| Asort_Carrera | 0.1373 | 0.1142 | 0.0750 | -0.0252 | 0.2684 |
| Asort_Estudios | 0.3134 | 0.1062 | 0.1580 | -0.1628 | 0.4385 |
| Asort_Genero | -0.0157 | -0.0418 | 0.1493 | -0.3498 | 0.2498 |
| Asort_Partido | 0.0096 | -0.0165 | 0.0782 | -0.1586 | 0.1432 |
| Densidad | 0.0466 | 0.0476 | 0.0069 | 0.0345 | 0.0623 |
| Dist_geodesica | 2.8486 | 3.2939 | 0.3193 | 2.8031 | 4.0272 |
| Transitividad | 0.6290 | 0.0833 | 0.0325 | 0.0228 | 0.1497 |
| Estadística | Observado | Media | Desv. Est. | IC Inferior | IC Superior |
|---|---|---|---|---|---|
| Asort_Carrera | 0.0471 | 0.0489 | 0.0367 | -0.0154 | 0.1292 |
| Asort_Estudios | 0.0983 | 0.0840 | 0.0888 | -0.0784 | 0.2676 |
| Asort_Genero | 0.2152 | 0.2063 | 0.0963 | 0.0351 | 0.3979 |
| Asort_Partido | 0.0506 | 0.0261 | 0.0422 | -0.0533 | 0.1133 |
| Densidad | 0.0505 | 0.0515 | 0.0050 | 0.0420 | 0.0614 |
| Dist_geodesica | 2.8897 | 2.9629 | 0.1472 | 2.7148 | 3.3038 |
| Transitividad | 0.6515 | 0.0918 | 0.0202 | 0.0514 | 0.1329 |
Por último, el modelo de sociabilidad con covariables nodales también muestra un buen ajuste en términos de los efectos homofílicos de las variables nodales incluidas en la formulación, así como un muy buen ajuste de la densidad de la red. EL comportamiento de la distancia geodésica promedio se repite ya que para las redes de Duque y Petro ésta es sobre-estimada por el modelo, mientras que es subestimada para la red de Santos.
Para todas las redes de pdoer, el modelo de sociabilidad no reproduce ni de cerca a la transitividad de los datos observados.
Conclusión: Haber agregado a las variables nodales dentro del modelo ERGM y el modelo de sociabilidad muestra una notable mejoría en la recuperación de los efectos homofílicos de estas covariables.
No obstante, no hay modelo perfecto ya que ninguno de estos logró capturar la alta transitividad de las redes de funcionarios. Probablemente, este desempeño hubiese sido mejor si en la formulación del modelo ERGM, por ejemplo, se hubiese incorporado un término para medir la transitividad.
Si tuviese que elegir uno de los tres modelos como el mejor para estas bases de datos, seleccionaría el ERGM.
Ahora, analizamos la red de poder de Colombia sin distinguir por gobierno. En esta base total contamos con 264 funcionarios públicos.
# Base completa
base_total <- poder %>%
mutate(EDAD = as.numeric(EDAD)) %>%
select(-GOBERNACION)
# Matriz bipartita
A_total <- crear_bipartita(base_total)
# Matriz de adyacencia binaria
Y_total <- (A_total %*% t(A_total) > 0) * 1
diag(Y_total) <- 0
dim(Y_total)## [1] 264 264
# RED
red_total <- graph_from_adjacency_matrix(Y_total, mode = "undirected", diag = FALSE)
V(red_total)$EDAD <- base_total$EDAD
V(red_total)$GENERO <- base_total$GENERO
V(red_total)$PARTIDO_POLITICO <- base_total$PARTIDO_POLITICO
V(red_total)$NIVEL_DE_ESTUDIOS <- base_total$NIVEL_DE_ESTUDIOS
V(red_total)$UNIVERSIDAD <- base_total$UNIVERSIDAD
V(red_total)$CARRERA <- base_total$CARRERAbrowsable(tagList(
tags$div(style = "margin-bottom:10px;", HTML("<b>Género:</b> ● Mujer ▲ Hombre")),
graficar_red_interactiva(red_total,'Red de poder en Colombia (2002-2026)', top_n = 0, centralidad = "eigenvector")
))## GERMAN VARGAS LLERAS MARTA LUCIA RAMIREZ BLANCO
## 0.4799270 0.4755877
## CARLOS HOLMES TRUJILLO GARCIA MAURICIO CARDENAS SANTAMARIA
## 0.4721724 0.4663121
## FEDERICO ALONSO RENJIFO VELEZ
## 0.4442568
## CARLOS HOLMES TRUJILLO GARCIA GERMAN VARGAS LLERAS
## 0.1695966 0.1069419
## MARTA LUCIA RAMIREZ BLANCO MAURICIO CARDENAS SANTAMARIA
## 0.1058832 0.1036709
## JUAN CARLOS ESGUERRA PORTOCARRERO
## 0.1000975
## GERMAN VARGAS LLERAS ARMANDO ALBERTO BENEDETTI VILLANEDA
## 1.0000000 0.9440143
## JUAN FERNANDO CRISTO BUSTOS LUIS FERNANDO VELASCO CHAVES
## 0.9277751 0.9277751
## NANCY PATRICIA GUTIERREZ CASTANEDA
## 0.9277751
German Vargas Lleras se corona como el funcionario más importante de la red de poder bajo dos de los tres criterios evaluados, además de ser el único que aparece en el top 5 para las tres centralidades. Su importancia se debe, entre otras cosas, a los cargos de alta relevancia que ocupó durante los gobiernos de Uribe y Santos. Marta Lucia Ramirez, Carlos Holmes Trujillo y Mauricio Cardenas también sobresalen en la red ya que están cercanos a varios otros nodos en la red y se encuentran entre varias relaciones de la misma. Curiosamente, Carlos Holmes y Mauricio Cardenas solo tuvieron cargos públicos durante 1 solo gobierno (Duque y Santos, respectivamente), lo que habla del impacto de sus gestiones en esos periodos para que resalten en la red completa.
En último, Armando Benedetti, Juan Fernando Cristo, Luis Fernando Velasco y Nacny Gutierrez son importantes en cuanto a que se roden de varios funcionarios igualmente centrales o importantes.
En el proximo grafico se resaltan los funcionarios más importantes mediante el criterio de intermediación.
## [1] 7
## [1] 3.02339
# Visualización
caminos <- distance_table(red_total)$res
names(caminos) <- 1:length(caminos)
barplot(
prop.table(caminos),
xlab = "Distancia geodésica",
ylab = "F. Relativa",
border = "grey",
col = "grey",
main = "Distribución de distancias geodésicas para la\n red de poder de Colombia (2002-2026)"
)El diametro de la red de poder toma un valor de 7, el cual podemos considerar como pequeño dada la cantidad de nodos presentes. Recordando que este valor se tenía para las redes filtradas por gobierno, el hecho de que se consideren todos los funcionarios juntos no aumenta las distancias entre los actores. En promedio, los funcionarios públicos se conectan mediante 3 o 4 intermediarios.
Adicionalmente, para la red general es más marcado que la distancia geodésica más frecuente es 3. Esto nos dice que la red no es extremadamente compacta, pero tampoco dispersa.
## [1] 23
## [1] 0.09
Nuevamente el clique máximo representa una porción pequeña del total de nodos, evidenciando que los vínculos institucionales se encuentran divididos en grupos aislados. El grupo más grande de funcionarios con cargos comunes es de 23.
## [1] 0.0580136
## [1] 0.6593249
# intransitividad local
hist(transitivity(red_total, type = "local"), xlab = "Transitividad local", ylab = "Frecuencia",
main = "Distribución de la transitividad local")La red de proyección presenta una densidad de 0.058, la cual es menor a la de las redes de Uribe y Santos, aunque mayor a las redes de Duque y Petro. Como se ha visto en análisis previos, esta no es una red que se caracterice por ser densa.
Por otro lado, en temas de transitividad se obtiene un valor de 0.659 el cual significa que aproximadamente el 65,9% de las triadas potencialmente cerrables están en efecto cerradas. En otras palabras, cuando dos funcionarios comparten vínculos con un mismo tercero, existe una alta probabilidad de que también se encuentren conectados entre sí. Esto sugiere la presencia de grupos cohesivos y patrones de interacción altamente estructurados. Este valor podría ser mayor, de no ser por la alta dispersión que se evidencia en las transitividades locales.
## [1] TRUE
## [1] 1
## [1] 1
A diferencia de las redes de los gobiernos, la red total si tiene la particularidad de ser conectada. Esta además tiene una conectividad nodal y d aristas de 1, lo que se traduce en que basta con remover una arista o un funcionariod e la red para que esta se desconecte
## [1] 1
##
## 264
## 1
## + 7/264 vertices, named, from a5e3778:
## [1] FRANCIA ELENA MARQUEZ MINA
## [2] GUILLERMO ANTONIO HERRERA CASTANO
## [3] CIELO ELAINNE RUSINQUE URREGO
## [4] MARIA DEL ROSARIO GUERRA DE LA ESPRIELLA
## [5] MAURICIO PERFETTI DEL CORRAL
## [6] AURORA VERGARA FIGUEROA
## [7] JUAN MANUEL SANTOS CALDERON
## [1] 0.02651515
Se cuentan con 7 nodos de articulación, los cuales representan el 2,65% de los actores totales. Algunos de los funcionarios que hacen parte de este grupo son Francia Marquez y Juan Manuel Santos.
# algoritmos
set.seed(123)
clust_fast_greedy <- cluster_fast_greedy(red_total)
clust_leading_eigen <- cluster_leading_eigen(red_total)
clust_walktrap <- cluster_walktrap(red_total)
clust_louvain <- cluster_louvain(red_total)
clust_label_prop <- cluster_label_prop(red_total)
#clust_optimal <- cluster_optimal(g2_bogota) no se icluye por costo computacional
clust_infomap <- cluster_infomap(red_total)tabla_modularidad <- data.frame(
Medida = c("Modularidad", "Número de clusters"),
FastGreedy = c(sprintf("%.3f", modularity(clust_fast_greedy)),
sprintf("%.0f", length(clust_fast_greedy)) ),
LeadingEigen = c(sprintf("%.3f", modularity(clust_leading_eigen)),
sprintf("%.0f", length(clust_leading_eigen)) ),
Walktrap = c(sprintf("%.3f", modularity(clust_walktrap)),
sprintf("%.0f", length(clust_walktrap)) ),
Louvain = c(sprintf("%.3f", modularity(clust_louvain)),
sprintf("%.0f", length(clust_louvain)) ),
LabelProp = c(sprintf("%.3f", modularity(clust_label_prop)),
sprintf("%.0f", length(clust_label_prop)) ),
Infomap = c(sprintf("%.3f", modularity(clust_infomap)),
sprintf("%.0f", length(clust_infomap)) ),
check.names = FALSE
)
kable(tabla_modularidad, align = "c",
caption = "Modularidad y número de clusters por algoritmo") %>%
kable_styling(bootstrap_options = c("bordered", "striped", "hover"),
full_width = FALSE)| Medida | FastGreedy | LeadingEigen | Walktrap | Louvain | LabelProp | Infomap |
|---|---|---|---|---|---|---|
| Modularidad | 0.632 | 0.628 | 0.642 | 0.671 | 0.631 | 0.644 |
| Número de clusters | 11 | 17 | 20 | 12 | 23 | 25 |
Al comparar 6 algoritmos de agrupamiento se obtiene que el mejor comportado es el algoritmo semi-jerárquico Louvain, con una modularidad de 0.671 y 12 clusters generados.
library(RColorBrewer)
# Comunidades de Louvain
V(red_total)$grupo <- factor(clust_louvain$membership)
# Nodos
nodes <- data.frame(
id = V(red_total)$name,
label = "",
grupo = V(red_total)$grupo,
partido = V(red_total)$PARTIDO_POLITICO,
edu = V(red_total)$NIVEL_DE_ESTUDIOS,
edad = V(red_total)$EDAD,
genero = V(red_total)$GENERO,
stringsAsFactors = FALSE
)
# Tamaño según edad
nodes$value <- rescale(nodes$edad, to = c(10, 40))
# Comunidad para colorear
nodes$group <- as.character(nodes$grupo)
# Forma según género
nodes$shape <- case_when(
nodes$genero == "F" ~ "dot",
nodes$genero == "M" ~ "triangle",
TRUE ~ "square"
)
# Información emergente al pasar el mouse
nodes$title <- paste0(
"<b>", nodes$id, "</b><br>",
"Partido: ", nodes$partido, "<br>",
"Nivel educativo: ", nodes$edu, "<br>",
"Edad: ", nodes$edad, "<br>"
)
# Aristas
edges <- igraph::as_data_frame(red_total, what = "edges")
names(edges)[1:2] <- c("from", "to")
edges$width <- 2
# Colores para las comunidades
comunidades <- sort(unique(nodes$group))
n_com <- length(comunidades)
if(n_com <= 8){
colores <- brewer.pal(max(3, n_com), "Set1")[1:n_com]
} else {
colores <- colorRampPalette(brewer.pal(8, "Set1"))(n_com)
}
# Objeto visNetwork
graf <- visNetwork(nodes,
edges,
width = "100%",
height = "850px",
main = "Agrupamiento para la Red de poder en Colombia")
# Colores y leyenda a cada comunidad
for(i in seq_along(comunidades)){
graf <- graf %>%
visGroups(groupname = comunidades[i], color = colores[i])
}
# Graficar
graf %>%
visEdges(smooth = FALSE,
color = list(color = "rgba(180,180,180,0.5)") ) %>%
visLegend(useGroups = TRUE,
position = "right",
main = "Comunidades") %>%
visOptions(highlightNearest = list(enabled = TRUE, hover = TRUE) ) %>%
visInteraction(dragNodes = TRUE,
dragView = TRUE,
zoomView = TRUE,
navigationButtons = TRUE) %>%
visPhysics(solver = "forceAtlas2Based",
stabilization = list(enabled = TRUE, iterations = 1000) ) %>%
visEvents(stabilizationIterationsDone ="function () {this.setOptions({physics:false});}" ) %>%
visLayout(randomSeed = 123)Visualmente el algoritmo asignó en los mismos clusters a funcionarios muy conectados entre si, es decir, hay evidentes grupos de funcionarios con coincidencias institucionales mutuas y el algoritmo está recuperando esta característica proveniente de la red bipartita (increible). Por ejemplo, al contrastar con la base de datos vemos que los funcionarios del cluster 12 han sido nombrados como directores del DANE (mismo cargo público). De esta forma, podemos hacernos una idea de quienes han ocupado que cargos desde 2002 hasta 2026 sin revisar individuo por individuo en la base.
vars_nominales <- c("GENERO", "PARTIDO_POLITICO", "NIVEL_DE_ESTUDIOS", "UNIVERSIDAD", "CARRERA")
vars_numericas <- c("EDAD")
tabla_asort <- data.frame(Red = "Total")
# Variables categóricas
for(var in vars_nominales){
tabla_asort[[var]] <- assortativity_nominal(red_total,
types = as.integer(factor(vertex_attr(red_total, var))),
directed = FALSE)
}
# Variables numéricas
for(var in vars_numericas){
g2 <- delete_vertices(red_total, V(red_total)[is.na(vertex_attr(red_total, var))])
tabla_asort[[var]] <- assortativity(g2,
values = vertex_attr(g2, var),
directed = FALSE)
}
kable(tabla_asort, digits = 3, align = "c",
caption = "Coeficientes de asortatividad") %>%
kable_styling(bootstrap_options = c("bordered"))| Red | GENERO | PARTIDO_POLITICO | NIVEL_DE_ESTUDIOS | UNIVERSIDAD | CARRERA | EDAD |
|---|---|---|---|---|---|---|
| Total | 0.122 | 0.049 | 0.095 | 0.03 | 0.098 | 0.137 |
Para la red total de poder del país, los nodos si tienden a relacionarse con otros que compartan sus características. En otras palabras, los funcionarios del país tienden a compartir cargos públicos con otros de edad similar, mismo género, mismo nivel de estudios, etc. La variable con mayor homofilia es la edad, seguida del género y la carrera estudiada.
El objetivo de este taller es analizar 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 anti
mafia “Montagna”. La investigación se centró en la familia “Mistretta” y
el clan “Batanesi”, que infiltraron diversas actividades económicas
entre 2003 y 2007. Los datos se encuentran en el archivo
mafia.RData. La red representa individuos y relaciones
asociadas con actividades ilícitas, donde los nodos corresponden a
personas o entidades y los arcos a relaciones de influencia o
interacción. La red es dirigida y ponderada.
Explorar la estructura de los datos de nodos y aristas, identificar la información relevante para el análisis, verificar si la red es simple y ponderada, construir la red, calcular su orden y tamaño, generar una visualización inicial e identificar los nodos más conectados.
Analizar la estructura global de la red mediante el grado nodal, la conectividad, los componentes, el diámetro y las distancias geodésicas, complementando el estudio con una visualización en la que los nodos se distingan por clan y tamaño según su grado.
Evaluar la importancia estructural de los nodos mediante grado, fuerza y medidas de centralidad, estudiar la cohesión de la red a través de cliques, densidad y transitividad, y apoyar la interpretación con visualizaciones adecuadas.
Aplicar y comparar métodos de detección de comunidades, evaluar la modularidad de las particiones, analizar la asortatividad por grado y por clan si corresponde, contrastar las comunidades detectadas con los clanes conocidos y visualizar la partición más adecuada.
Ajustar los modelos considerados previamente incorporando la información nodal relevante cuando sea pertinente, presentar e interpretar sus resultados mediante tablas y gráficos, evaluar su bondad de ajuste con base en densidad, transitividad, asortatividad y distancia geodésica promedio, y comparar su capacidad predictiva mediante validación cruzada de 5 particiones. Las secciones de bondad de ajuste y capacidad predictiva deben tener una extensión máxima de 1000 palabras cada una.
Solución Punto 3
## [1] 143
## [1] 325
## [1] TRUE
## [1] TRUE
Esta es una red no dirigida, simple y ponderada. Cuenta con 143 actores y 325 aristas entre ellos.
V(G)$grado <- degree(G)
set.seed(123)
g_tbl <- as_tbl_graph(G)
layout <- create_layout(g_tbl, layout = "igraph", algorithm = "nicely")
ggraph(layout) +
geom_edge_link(color = "grey60",
alpha = 0.4,
aes(width = weight),
end_cap = circle(0, "mm") ) +
geom_node_point(color = "darkblue",
aes(size = grado),
alpha = 0.7
) +
scale_edge_width(range = c(1, 6), guide = "none") +
scale_size(range = c(2,8), guide = "none") +
theme_void() +
theme(legend.position = "right") +
labs(title = "Red de mafia"#,
#subtitle = "El tamaño de los vértices es proporcional a la edad"
)## N18 N47 N68 N27 N61
## 40 26 25 19 18
Los actores más conectados de la red son aquellos que tienen un mayor grado. El top 5 de nodos más conectados son N18, N47, N68 y N61. Para entenderlos mejor basta con mirar sus roles. N18 hace parte de la cúpula ejecutiva de la familia Mistretta, N47 es el jefe adjunto de la familia Batanesi, N68 y N27 son parte de la cúpula ejecutiva de esta, y finalmente N61 hace parte de la cúpula de la familia Mistretta.
Vemos que los actores con mayor grado tienen posiciones altas dentro de las dos familias en las que se centró la investigación.
## N47 N18 N68 N27 N29
## 100 86 55 48 43
Cuando el criterio ahora es la fuerza, vemos que el top 5 varía levemente. En primera posición con 100 interacciones ilícitas tenemos a N47, el jefe adjunto de la familia Batanesi. Los demás particapntes de este top hacen parte de las dos familias ya mencionadas: Batanesi y Mistretta, además de un emprendedor/comerciante identificado como el nodo N29.
A continuación se lleva a cabo un análisis local y estructural de la red de mafia. Primero, generamos una visualización considerando los roles de los individuos, en cuanto a las familias a las que pertenecen. Individuo que en su rol no tenga asociada una de las familias, se le asigna un NA.
library(stringr)
familias <- c("Batanesi", "Mistretta", "Mazzaroti", "Caltagirone", "Tortorici", "enterpreneur")
V(G)$familia <- str_extract(V(G)$role, paste(familias, collapse = "|") )
V(G)$familia[V(G)$familia == ""] <- NA
# Para las aristas
familias_mafia <- c("Batanesi", "Mistretta", "Mazzaroti", "Caltagirone", "Tortorici")
V(G)$familia_mafia <- ifelse(V(G)$familia %in% familias_mafia, V(G)$familia, NA)
ed <- ends(G, E(G), names = FALSE)
fam1 <- V(G)$familia_mafia[ed[,1]]
fam2 <- V(G)$familia_mafia[ed[,2]]
E(G)$tipo_arista <- NA
E(G)$tipo_arista[!is.na(fam1) & !is.na(fam2) & fam1 == fam2] <- "Misma familia"
E(G)$tipo_arista[!is.na(fam1) & !is.na(fam2) & fam1 != fam2] <- "Entre familias"
E(G)$tipo_arista[xor(!is.na(fam1), !is.na(fam2))] <- "Con intermediarios"
E(G)$tipo_arista[is.na(fam1) & is.na(fam2)] <- "Sin familia"
set.seed(123)
g_tbl <- as_tbl_graph(G)
layout <- create_layout(g_tbl, layout = "igraph", algorithm = "nicely")
ggraph(layout) +
geom_edge_link(aes(width = weight),
color = "grey60",
alpha = 0.4,
end_cap = circle(0, "mm") ) +
geom_node_point(aes(color = familia,
size = grado),
alpha = 0.7
) +
scale_edge_width(range = c(1, 6), guide = "none") +
scale_size(range = c(2,8), guide = "none") +
scale_color_brewer(palette = "Set1", name = "Familia/Clan", na.value = "grey50") +
theme_void() +
theme(legend.position = "right") +
labs(title = "Red de mafia",
subtitle = "El tamaño de los vértices es proporcional al grado"
)Note que, como vimos en el análisis previo de los grados y las fuerzas, los nodos más grandes (i.e. con mayores interacciones ilegales) hacen parte del clan Batanesi y la familia Mistretta. Esta última parece haber tenido muchas relaciones con varios individuos que no necesariamente hacían parte de su misma familia, mientras que los Bataseni están muy conectados entre ellos de forma que se agrupan en una parte específica del grafo.
Aparte de las familias, los emprendedores (enterpreneur) también juegan un papel importante de puentes entre varios nodos.
set.seed(123)
g_tbl <- as_tbl_graph(G)
layout <- create_layout(g_tbl, layout = "igraph", algorithm = "nicely")
ggraph(layout) +
geom_edge_link(aes(width = weight,
color = tipo_arista),
alpha = 0.4,
end_cap = circle(0, "mm") ) +
geom_node_point(aes(size = grado),
color = "grey40",
alpha = 0.5
) +
scale_edge_width(range = c(1, 6), guide = "none") +
scale_size(range = c(2,8), guide = "none") +
scale_edge_color_manual(values = c("Entre familias" = "red",
"Con intermediarios" = "orange",
"Misma familia" = "darkgreen",
"Sin familia" = "grey70"),
name = "Tipo de vínculo") +
theme_void() +
theme(legend.position = "right") +
labs(title = "Red de mafia",
subtitle = "El tamaño de los vértices es proporcional al grado"
)##
## Con intermediarios Entre familias Misma familia Sin familia
## 157 19 57 92
Dada la clasificación natural de los nodos en sus clanes/familias, también se puede pensar en clasificar las aristas dependiendo del tipo de nodos que se conectan. En este caso, se decide explorar los tipos de enlaces que se dan dentro de la misma familia, entre familias, entre una familia y un intermediario, y finalmente, las que son entre personas sin rol en alguna familia.
Note primero que las interacciones entre miembros de la mafia e intermediarios (no miembros de estas familias) corresponde al 48.31% de las conexiones totales de la red. Dejando ver la importancia y concurrencia de estos individuos en una red de actividades ilícitas, ya que juegan el papel de puentes entre las interacciones orquestadas por los grandes cabecillas. Los vínculos entre familias son menos frecuentes, dejando ver que son escasos los intercambios ilícitos directos entre miembros de grandes carteles.
## [1] 12
## [1] 4.178872
# Visualización
caminos <- distance_table(G)$res
names(caminos) <- 1:length(caminos)
barplot(
prop.table(caminos),
xlab = "Distancia geodésica",
ylab = "F. Relativa",
border = "grey",
col = "grey",
main = "Distribución de distancias geodésicas"
)La distancia geodésica entre dos vértices es la longitud del camino más corto que los conecta. Un diámetro de red igual a 12, como valor máximo de estas distancias, nos dice que los dos individuos más alejados necesitan 12 interacciones delictivas para conectarse por el camino más corto.
Por su parte, la distancia geodésica promedio nos dice que se requieren de 4 a 5 personajes intermedios en promedio para conectar a dos actores. Complementariamente, la distribución de las distancias geodésicas muestra un comportamuiento acampanado agrupado alrededor de 3.
Se calcula el top 5 de individuos más centrales según los criterios vistos en clase.
Nota: La centralidad en términos de grado y fuerza ya se trató previamente.
## N16 N73 N106 N2 N17
## 1 1 1 1 1
## N18 N47 N68 N61 N27
## 0.4654077 0.1991987 0.1944497 0.1558615 0.1190395
## N47 N45 N27 N51 N48
## 1.0000000 0.5910204 0.5122061 0.5044318 0.4898131
Los actores más importantes de la red de mafia bajo el criterio de cercanía son, curiosamente, todos personajes no identificados dentro de algún clan/familia, personas sin rol. Esto nos deja ver que aunque en movimientos ilícitos los grandes capos son la base, la red se mueve con ayuda de varios actores que no necesariamente tienen el mismo ‘prestigio’ o origen.
La centralidad por intermediación nos devela un comportamiento diferente al anterior, ya que acá si priman los miembros de los grandes clanes/familias. Encabezando el top, está un miembro de la familia Mistretta (N18) seguido de 2 miembros del clan Batanesi (N47 y N68). Este hallazgo complementa casi a la perfección lo mencionado anteriormente: los actores más cercanos en una red de mafia pueden ser personas ‘sin mayor importancia’, pero los grandes carteles siempre estarán entre esas conexiones ilícitas.
Por último, miembros del clan Bataseni son los más importantes de la red en cuanto a que se rodean de varios actores igualmente importantes. Evidenciando que este clan tiene una preferencia por interactuar con actores igualmente conocidos en el “negocio”.
# frecuencias de clanes
table(sapply(X = cliques(graph = G), FUN = length)) # de tamaño 7, hay 1 clique##
## 1 2 3 4 5 6 7
## 143 325 315 191 66 12 1
## [[1]]
## + 7/143 vertices, named, from eb7f679:
## [1] N8 N3 N9 N7 N6 N5 N4
## [1] 7
En la red de mafia contamos con clanes de siete tamaños diferentes. El clan/clique más grande se comprende de solo 7 individuos. Este clan máximo tiene a un miembro de la familia Caltagirone, un emprendedor y demás actores sin rol específico. Que el númeor clan tome el valor de 7 nos indica que las interacciones ilegales suelen darse en grupos completamente conectados pequeños.
## [1] 0.03201024
# Puedo hacer la densidad para los subgrafos inducidos por Batanesi y Mistretta
ids_batanesi <- which(V(G)$familia == "Batanesi")
vecinos <- unique(unlist(neighborhood(G, order = 1, nodes = ids_batanesi)))
g_batanesi <- induced_subgraph(G, vecinos)
ids_mistre <- which(V(G)$familia == "Mistretta")
vecinos <- unique(unlist(neighborhood(G, order = 1, nodes = ids_mistre)))
g_mistre <- induced_subgraph(G, vecinos)
edge_density(graph = g_batanesi)## [1] 0.1062194
## [1] 0.05809802
Un 0.032 muestra que esta es una red poco densa, bastante desconectada y considerablemente lejana de volverse un clan.
Se decide mostrar adicionalmente la densidad de las redes inducidas por los vecindarios de los miembros de los dos grupos más importantes ya mencionados, Batanesi y Mistretta. Observe que el contexto relacional del clan Bitanesi es más denso que la red completa y que la red de vecinos de la familia Mistretta.
## [1] 0.2785146
# intransitividad local
hist(transitivity(G, type = "local"), xlab = "Transitividad local", ylab = "Frecuencia",
main = "Distribución de la transitividad local")Estamos ante una red moderadamente transitiva. Era de esperarse que si un individuo A tiene una conexión ilícita con B, y B también tiene una conexión de este tipo con C, no es común que A y C vayan a tener una conexión de este tipo.
La distribución de las transitividades locales muestra una gran dispersión en esta variable. Hay 40 relaciones 100% transitivas, mientras que las demás tienen un comportamiento casi uniforme en el intervalo [0, 1], sin superar una frecuencia de 10 cada una. Luego de la transitividad de 1 la más frecuente es 0.5.
## [1] FALSE
## [1] 0
## [1] 0
Como habíamos visto en las visualizaciones, esta red no está conectada, por lo que sus conectividades nodal y de aristas toman el valor de cero.
Para estudiar mejor la red, procedemos a encontrar la componente gigante de esta.
## [1] 5
##
## 2 3 134
## 3 1 1
La red de mafia cuenta con 5 componentes. Una de ellas es una tríada, hay otras tres aisladas que son simplemente díadas; y una última que se compone de 134 nodos. Esta última es, en evidencia, la componente gigante.
## [1] 1
## [1] 1
## + 15/134 vertices, named, from 77de329:
## [1] N14 N77 N75 N109 N99 N23 N18 N56 N52 N36 N61 N47 N68 N27 N11
## [1] 0.1119403
La componente gigante tiene conectividad nodal y de aristas de 1. Es decir que basta con remover un nodo o una arista de esta red, para que se desconecte. Hay 15 actores que se consideran puntos de articulación, los cuales corresponden al 11,19% del orden de la red.
Se evalua la homofilia de los nodos en cuanto a su grado y al rol de estos en la mafia.
## [1] -0.2325615
Los coeficientes de asortatividad miden el nivel de homofilia del grafo a partir de alguna característica de los vértices.Que la asortatividad del grado sea negativa deja ver una falta de homofilia en esta característica, es decir, que los nodos tienden a relacionarse con otros que tienen un número de conexiones distinto.- En el contexto de nuestra red esto se traduce a que los grandes mafiosos tienden a relacionarse ilegalmente con personajes “pequeños”.
V(G)$familia <- ifelse(is.na(V(G)$familia), "Sin_familia", V(G)$familia)
clan <- as.numeric(as.factor(V(G)$familia))
assortativity_nominal(G, types = clan)## [1] 0.1573403
Del mismo modo, se puede cuantificar la homofilia inducida por la característica nodal del rol. En este caso, el coeficiente de asortatividad toma un valor positivo, lo que quiere decir que los individuos con el mismo rol tienden a interactuar entre si, donde el rol es la familia a la que hacen parte (recordando que si no pertenecen a ninguna familia, esta variable se toma como ‘Sin familia’).
Se corren 6 algoritmos de agrupamiento sobre la red de mafia.
set.seed(123)
# algoritmos
clust_fast_greedy <- cluster_fast_greedy(G)
clust_leading_eigen <- cluster_leading_eigen(G)
clust_walktrap <- cluster_walktrap(G)
clust_louvain <- cluster_louvain(G)
clust_label_prop <- cluster_label_prop(G)
#clust_optimal <- cluster_optimal(G)
clust_infomap <- cluster_infomap(G)# gráficos
igraph_options(vertex.size = 10, vertex.frame.color = "black")
par(mfrow = c(3, 3), mar = c(0, 0, 2, 0))
# extraer coordenadas del layout de ggraph
layout_mat <- as.matrix(layout[, c("x", "y")])
plot(G, vertex.label = NA, layout = layout_mat, vertex.color = clust_fast_greedy$membership,
main = paste0("fast greedy: ", "Mod = ", round(modularity(clust_fast_greedy), 4)))
plot(G, vertex.label = NA, layout = layout_mat, vertex.color = clust_leading_eigen$membership,
main = paste0("leading eigen: ", "Mod = ", round(modularity(clust_leading_eigen), 4)))
plot(G, vertex.label = NA, layout = layout_mat, vertex.color = clust_walktrap$membership,
main = paste0("walktrap: ", "Mod = ", round(modularity(clust_walktrap), 4)))
plot(G, vertex.label = NA, layout = layout_mat, vertex.color = clust_louvain$membership,
main = paste0("louvain: ", "Mod = ", round(modularity(clust_louvain), 4)))
plot(G, vertex.label = NA, layout = layout_mat, vertex.color = clust_label_prop$membership,
main = paste0("label prop: ", "Mod = ", round(modularity(clust_label_prop), 4)))
plot(G, vertex.label = NA, layout = layout_mat, vertex.color = clust_infomap$membership,
main = paste0("infomap: ", "Mod = ", round(modularity(clust_infomap), 4)))De los algoritmos de agrupamiento estudiados, el que maximiza la modularidad es Fast Greedy, un algoritmo jerárquico que directamente busca maximizar esta estadística.
Es por esto que a continuación se estudian sus comunidades en el siguiente gráfico.
# agregar comunidades como atributo de nodo
g_tbl <- g_tbl %>%
mutate(comunidad = factor(clust_fast_greedy$membership))
set_graph_style(plot_margin = margin(1,1,1,1))
set.seed(123)
layout <- create_layout(g_tbl, layout = "igraph", algorithm = "nicely")
ggraph(layout) +
geom_edge_link(aes(width = weight),
alpha = 0.4,
color = "gray60",
end_cap = circle(0, "mm") ) +
geom_node_point(aes(color = comunidad,
size = grado),
alpha = 0.7) +
scale_edge_width(range = c(1, 6), guide = "none") +
scale_size(range = c(2, 8), guide = "none") +
theme_void() +
labs(
title = "Agrupamiento de la red de mafia mediante el algoritmo Fast Greedy",
subtitle = "Modularidad = 0.5717",
color = "Comunidad"
)EL algoritmo devuelve 10 comunidades que parecen repartirse de acuerdo a las vecindades de los nodos más conectados de la red, ya que observamos al menos un nodo de alto grado hace parte de cada comunidad. Como es comun, las 4 componentes aisladas de la componente gigante son una comunidad cada una.
Recordando la clasificación de los individuos de acuerdo a sus familias/clanes vemos que la comunidad dos tiene fuertes similitudes con el grupo de miembros del clan Batanesi.
## [1] TRUE
## Community sizes
## 1 2 3 4 5 6 7 8 9 10
## 33 35 23 19 7 17 3 2 2 2
El cluster más grande se comprende de 35 actores, el cual como ya se mencionó, abarca a los miembros del clan Batanesi. El que le sigue en tamaño es el cluster 1 con 33 individuos, el cual se ubica en la parte superior izquierda de la componente gigante y se caracteriza por contener al nodo más conectado de la red.
Que hayan tantos clusters nos habla de una heterogeneidad de los datos, aspecto ya estudiado en los análisis previos.
Finalmente, dada la agrupación natural de los nodos en sus roles en la mafia, se compara con esta la agrupación generada por el algoritmo mediante el indice de rand ajustado
Vease que comparando la agrupación del algoritmo con la inducida por el rol de los individuos se obtiene un resultado interesante. Por un lado, el indice de rand dice que el 59.5% de las díadas reciben la misma clasificación en ambas particiones. No obstante, en redes con varias comunidades o grupos muy desvalanceados, como la nuestra, gran parte de esa coincidencia viene del azar. El indice de rand ajustado, al corregir este error, demuestra que las comunidades encontradas por Fast Greedy son apenas mejores que una partición aleatoria para recuperar las divisiones de las familias mafiosas.
Antes de ajustar los modelos, veamos la matriz de adyacencia de la red observada.
Nota: Recuérdese que en estos modelos se necesita una red binaria.
#heat.plot0(as.matrix(as_adjacency_matrix(G)), labs = NA)
Y <- as.matrix(as_adjacency_matrix(graph = G, names = F))
corrplot::corrplot(corr = Y,
col.lim = c(0,1),
method = "color",
tl.pos = "n",
cl.pos = "r",
mar = c(0, 0, 3, 0),
col = colorRampPalette(c("blue","white","red"))(200))
title(main = "Matriz de adyacencia")En el modelo de Erdös-Rényi todas las probabilidades de interacción son las mismas: \(\theta\). Dado que este es en modelo muy sencillo, se sabe que en este caso la estimación por máxima verosimilitud está dada por \(\hat{\theta}_{MLE}=\text{den}(\mathbf{Y})\).
## [1] 0.03201024
# log-verosimilitud
n <- vcount(G)
m <- (n*(n-1))/2
s <- m*edge_density(G, loops = FALSE)
loglik <- function(theta) s*log(theta) + (m-s)*log(1-theta)
# gráfico
par(mfrow = c(1,2), mar = c(3,3,1,1), oma = c(0,0,3,0), mgp = c(1.75,0.75,0) )
curve(expr = loglik(x), from = 0, to = 0.5, lwd = 2, xlab = expression(theta), ylab = "Log-verosimilitud")
abline(v = theta_hat, col = 2, lty = 2)
curve(expr = loglik(x), from = 0.02, to = 0.04, lwd = 2, xlab = expression(theta), ylab = "Log-verosimilitud")
abline(v = theta_hat, col = 2, lty = 2)
mtext("Función de log-verosimilitud del parámetro", side = 3, outer = TRUE, line = 1, cex = 1.2)Este modelo agrega que las redes cumplan con una secuencia de grados dada.
## [1] 143
## [1] 325
## N0 N1 N3 N11 N12 N4 N5 N28 N25 N6 N7 N8 N10 N13 N14 N43
## 1 1 8 15 16 6 10 5 13 10 6 6 5 6 7 9
## N27 N29 N64 N47 N68 N45 N54 N36 N48 N16 N18 N100 N19 N22 N23 N26
## 19 16 7 26 25 14 7 10 12 1 40 4 11 17 9 2
## N32 N33 N34 N61 N70 N39 N40 N41 N46 N49 N50 N51 N52 N56 N63 N73
## 3 4 6 18 7 7 4 3 4 5 8 11 4 4 5 1
## N75 N76 N79 N80 N81 N84 N85 N86 N89 N91 N93 N101 N102 N103 N104 N108
## 11 7 5 4 3 4 8 4 11 2 5 1 1 1 1 1
## N127 N132 N125 N140 N146 N153 N110 N109 N105 N117 N128 N129 N130 N131 N133 N134
## 1 1 2 1 1 1 2 3 1 1 1 1 1 1 1 1
## N106 N112 N113 N114 N115 N116 N118 N119 N126 N136 N138 N121 N123 N147 N151 N139
## 1 1 1 2 1 1 1 1 1 1 1 1 1 1 2 2
## N149 N150 N152 N148 N124 N69 N137 N143 N2 N9 N21 N87 N83 N97 N17 N20
## 1 1 1 1 1 4 1 2 2 6 2 4 2 4 1 2
## N31 N35 N59 N60 N65 N71 N88 N30 N72 N77 N78 N98 N99 N24 N92 N62
## 2 5 2 2 1 2 2 2 2 6 2 3 4 1 2 1
## N66 N67 N42 N55 N90 N94 N95 N96 N57 N74 N82 N107 N111 N144 N145
## 3 1 3 2 2 2 4 3 1 1 3 1 1 2 1
## [1] 4.545455
La red de mafia cuenta con una secuencia de grados de longitud 143, cuya media es de 4.54 conexiones.
Ahora, vamos a obtener las probabilidades de interacción de este modelo.
# probabilidades de interacción
# Modelo 1: Modelo de Erdos-Renyi
# Modelo 2: Modelo de grafos aleatorios generalizado
set.seed(123)
B <- 1000
IP1 <- matrix(data = 0, nrow = n, ncol = n)
IP2 <- matrix(data = 0, nrow = n, ncol = n)
for (b in 1:B) {
IP1 <- IP1 + as.matrix(as_adjacency_matrix(sample_gnm(n = n, m = m, directed = F, loops = F)))/B
IP2 <- IP2 + as.matrix(as_adjacency_matrix(sample_degseq(out.deg = grados, method = "vl")))/B
}
# Visuazalición
par(mfrow = c(1,2), mar = c(1,1,1,1), oma = c(0,0,1,0), mgp = c(1.75,0.75,0) )
corrplot::corrplot(corr = IP1,
col.lim = c(0,1),
method = "color",
tl.pos = "n",
cl.pos = "r",
col = colorRampPalette(c("blue","white","red"))(200))
title(main = "M. Erdos-Renyi", line = -1)
corrplot::corrplot(corr = IP2,
col.lim = c(0,1),
method = "color",
tl.pos = "n",
cl.pos = "r")
title(main = "M. Grafos Aleatorio Generalizado", line = -1)Vea que el modelo de grafos aleatorios genera probabilidades de interacción mucho más homogéneas que varían muy poco del valor de la densidad de la red. Por su parte, la matriz de adyacencia obtenida por el modelo de grafos aleatorios generalizado trata de emular la alta conectividad de la matriz original en su parte superior izquierda. Naturalmente, observamos que las probabilidades de interacción ya no son constantes.
El modelo de grafos aleatorios exponencial estaba pensado para correrse inicialmente teniendo en cuenta la formación de triángulos, no obstante al hacer varias pruebas mostró problemas de convergencia por lo que al final nos quedamos con el número de aristas y el fecto de la covariable ‘familia’.
grafo <- asNetwork(G)
ergm_model<- formula(grafo ~ edges + nodematch("familia", diff = TRUE)) # + gwesp(1, fixed = TRUE)
summary(ergm_model)## edges nodematch.familia.Batanesi
## 325 46
## nodematch.familia.Caltagirone nodematch.familia.enterpreneur
## 0 18
## nodematch.familia.Mazzaroti nodematch.familia.Mistretta
## 1 10
## nodematch.familia.Sin_familia nodematch.familia.Tortorici
## 43 0
## Observed statistic(s) nodematch.familia.Caltagirone and nodematch.familia.Tortorici are at their smallest attainable values. Their coefficients will be fixed at -Inf.
## Starting maximum pseudolikelihood estimation (MPLE):
## Obtaining the responsible dyads.
## Evaluating the predictor and response matrix.
## Maximizing the pseudolikelihood.
## Finished MPLE.
## Evaluating log-likelihood at the estimate.
## Call:
## ergm(formula = ergm_model)
##
## Maximum Likelihood Results:
##
## Estimate Std. Error MCMC % z value Pr(>|z|)
## edges -3.40200 0.07065 0 -48.151 <1e-04 ***
## nodematch.familia.Batanesi 2.73083 0.19453 0 14.038 <1e-04 ***
## nodematch.familia.Caltagirone -Inf 0.00000 0 -Inf <1e-04 ***
## nodematch.familia.enterpreneur 0.18982 0.25057 0 0.758 0.4487
## nodematch.familia.Mazzaroti 2.70886 1.22678 0 2.208 0.0272 *
## nodematch.familia.Mistretta 2.14924 0.36546 0 5.881 <1e-04 ***
## nodematch.familia.Sin_familia -0.85575 0.16905 0 -5.062 <1e-04 ***
## nodematch.familia.Tortorici -Inf 0.00000 0 -Inf <1e-04 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Null Deviance: 14075 on 10153 degrees of freedom
## Residual Deviance: 2660 on 10145 degrees of freedom
##
## AIC: 2672 BIC: 2715 (Smaller is better. MC Std. Err. = 0)
## Warnings:
##
## * The following terms have infinite coefficient estimates due to an
## extreme sufficient statistic:
##
## nodematch.familia.Caltagirone, nodematch.familia.Tortorici
El coeficiente estimado del término edges de -3.402
indica que si dos actores no hacen parte de la misma familia/clan de
mafia tienen una probabilidad basal de conexión del 3.22%, la cual
podemos considerar como baja.
Sabiendo que el clan Batanesi y la familia Mistretta son de interés, analizamos sus coeficientes estimados. Si dos individuos son miembros de Batanesi, la probabilidad de que tengan conexiones ilícitas se posiciona en un 33.82%. Por su parte, si dos actores son miembros de Mistretta la probabilidad de conexión también aumenta y toma el valor de 22.22%.
Esto nos muestra que la pertenencia a estas grandes familias/clanes tiene un impacto positivo en la formación de aristas.
Bajo este modelo, la propensión de formar aristas entre nodos está influenciada por un agrupamiento de la red en distintas clases.
set.seed(123)
modelo_sbm <- blockmodels::BM_bernoulli(membership_type = "SBM_sym",
adj = Y,
verbosity = 0,
plotting = "")
# estimación
modelo_sbm$estimate()## [1] 4
plot(modelo_sbm$ICL, type = "b", pch = 16, xlab = "Q", ylab = "ICL", main = "Número óptimo de grupos")
abline(v = Q, col = "red", lty = 2)# Probabilidades de pertenencia
Z <- modelo_sbm$memberships[[Q]]$Z
# Asignaciones a los clusters
labs <- apply(Z, 1, which.max)
table(labs)## labs
## 1 2 3 4
## 21 93 23 6
## labs
## 1 2 3 4
## 0.14685315 0.65034965 0.16083916 0.04195804
El SBM muestra que el número óptimo de clases para la red de mafia es de 4. Estos cuatro grupos no están divididos homogéneamente en cuanto al tamaño, ya que vemos un grupo dominante con 93 nodos, y al mismo tiempo otro que consta de tan solo 6 individuos. Con estos valores, obtenemos los tamaños relativos de clases, los cuales son de 14,6%, 65%, 16,1% y 4,2% aproximadamente.
## [,1] [,2] [,3] [,4]
## [1,] 0.2175421131 0.015218094 0.0004895112 0.3556367
## [2,] 0.0152180941 0.006036326 0.0061734653 0.1165889
## [3,] 0.0004895112 0.006173465 0.2825078013 0.1221838
## [4,] 0.3556366661 0.116588877 0.1221838349 0.5196543
corrplot::corrplot(Pi, type = "full", col.lim = c(0,1), method = "shade",
addgrid.col = "gray90", tl.col = "black", mar = c(0,0,2,0),
title = "Probabilidades de interacción")Con este modelo también podemos analizar las probabilidades de interacción entre nodos de las distintas clases consideradas. Vemos que los nodos del grupo 4 tienen una alta probabilidad de interactuar entre si (52%), así como una probabilidad moderada-alta de interactuar con el grupo 1 (35%). Curiosamente, el modelo nos muestra que los nodos del grupo 2 casi no interactuan entre si.
# Ordenar la matriz de adyacencia
tmp <- get_adjacency_ordered(xi = labs,
A = Y)
cols <- RColorBrewer::brewer.pal(length(unique(labs)), "Set1")
par(mfrow = c(1,2), mar = c(1,1,1,1), oma = c(0,0,3,0))
set.seed(123)
plot(G, layout = layout_nicely,
vertex.label = NA,
vertex.size = 5,
vertex.color = cols[labs],
vertex.frame.color = cols[labs],
edge.color = adjustcolor("black", 0.1))
heat.plot0(mat = tmp$A, tick = FALSE, labs = NA, asp = 1)
abline(v = tmp$d + 0.5,
h = tmp$d + 0.5)
mtext("Agrupamiento de la red de mafia - Modelo de Bloques Estocásticos",
side = 3, outer = TRUE, line = 1, cex = 1.2)La visualización del agrupamiento de la red mediante SBM muestra características interesantes. Primero, confirmando los análisis previos, vemos que la matriz de adyacencia ordenada nos confirma que los nodos del grupo 2 están apenas conectados, cobijando a las componentes aisladas de la red y a los actores más periféricos y poco conectados de la componente gigante. El grupo 1 (de color rojo en la red) se ubica en el centro de la red repartido entre varias conexiones importantes de esta.
Resalta de este agrupamiento el grupo 4 (de color morado en el gráfico). A simple vista no es la gran cosa al contar solo con 6 individuos; no obstante, estos individuos son de los más conectados en la red, y por ende de los más importantes en el sistema de interacciones ilegales entre mafias.
Por último, las comunidades formadas por este modelo no tienen mucha similitud con las inducidas por el Rol de los personajes.
El modelo de distancia está construido bajo la idea de homofilia, ya que si dos individuos \(i\) y \(j\) están cerca en el espacio latente, su \(\theta_{ij}\) se mantiene.
En este modelo se van a incorporar dos dimensiones latentes.
library(latentnet)
modelo_ld <- ergmm(formula = grafo ~ euclidean(d = 2, G = 0),
seed = 123,
control = ergmm.control(burnin = 1000, # Iteraciones de burn-in
interval = 10, # Thinning: guardar 1 muestra cada 5 iteraciones
sample.size = 5000 # Total de muestras efectivas para inferencia
)
)
summary(modelo_ld)##
## ==========================
## Summary of model fit
## ==========================
##
## Formula: grafo ~ euclidean(d = 2, G = 0)
## Attribute: edges
## Model: Bernoulli
## MCMC sample of size 5000, draws are 10 iterations apart, after burnin of 1000 iterations.
## Covariate coefficients posterior means:
## Estimate 2.5% 97.5% 2*min(Pr(>0),Pr(<0))
## (Intercept) 0.206053 -0.068055 0.4962 0.144
##
## Overall BIC: 3096.747
## Likelihood BIC: 1959.643
## Latent space/clustering BIC: 1137.104
##
## Covariate coefficients MKL:
## Estimate
## (Intercept) -1.351423
## [1] 0.2060534
## [1] 0.5513319
El modelo latente de distancia no es significativo ya que su único coeficiente estimado (el intercepto) tiene un valor-p mayor al nivel de significancia. En cuanto al coeficiente estimado, que este tome el valor de 0.206 se traduce en una probabilidad basal de interacción en escala real del 55,1%.
Continuando con el análisis del modelo, evaluamos la convergencia de este.
# Datos simulados
x <- c(modelo_ld$sample$lpY)
iterations <- 1:length(x)
# Crear un dataframe para ggplot
data <- data.frame(Iteration = iterations, LogLikelihood = x)
# Calcular estadísticas
mean_x <- mean(x)
quantiles_x <- quantile(x, c(0.025, 0.975))
# Generar el gráfico
ggplot(data, aes(x = Iteration, y = LogLikelihood)) +
geom_point(alpha = 0.3, size = 0.5, color = "black") + # Puntos con transparencia
geom_hline(yintercept = mean_x, color = "blue", linetype = "dashed", linewidth = 1) + # Línea de la media
geom_hline(yintercept = quantiles_x, color = "red", linetype = "dotted", linewidth = 1) + # Líneas de los cuantiles
labs(x = "Iteración", y = "Log-verosimilitud", title = "") +
theme_minimal() +
theme(plot.title = element_text(hjust = 0.5)) # Centrar el título si es necesarioVemos que la cadena de la Log-verosimilitud se ve estacionaria alrededor de una constante, por lo que se afirma que el modelo converge correctamente.
# Datos simulados (reemplaza con fit$sample$beta)
x <- c(modelo_ld$sample$beta)
iterations <- 1:length(x)
data_chain <- data.frame(Iteration = iterations, Beta = x)
# Estadísticas
mean_x <- mean(x)
quantiles_x <- quantile(x, c(0.025, 0.975))
# Gráfico de la cadena
p_chain <- ggplot(data_chain, aes(x = Iteration, y = Beta)) +
geom_point(alpha = 0.3, size = 0.5, color = "black") + # Puntos con transparencia
geom_hline(yintercept = mean_x, color = "blue", linetype = "dashed", linewidth = 1) + # Línea de la media
geom_hline(yintercept = quantiles_x, color = "red", linetype = "dotted", linewidth = 1) + # Líneas de los cuantiles
labs(x = "Iteración", y = expression(beta), title = "Cadena") +
theme_minimal()
# Histograma de la distribución marginal
data_hist <- data.frame(Beta = x)
p_hist <- ggplot(data_hist, aes(x = Beta)) +
geom_histogram(aes(y = ..density..), bins = 30, fill = "gray90", color = "gray90") + # Histograma
geom_vline(xintercept = mean_x, color = "blue", linetype = "dashed", linewidth = 1) + # Línea de la media
geom_vline(xintercept = quantiles_x, color = "red", linetype = "dotted", linewidth = 1) + # Líneas de los cuantiles
labs(x = expression(beta), y = "Densidad", title = "Distr. marginal") +
theme_minimal()
library(grid)
# Combinar los gráficos
grid.arrange(p_chain, p_hist, ncol = 2,
top = textGrob("Inferencia sobre el intercepto",
gp = gpar(fontsize = 16, fontface = "plain") )
)## Warning: The dot-dot notation (`..density..`) was deprecated in ggplot2 3.4.0.
## ℹ Please use `after_stat(density)` instead.
## This warning is displayed once per session.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
La gráfica de la distribución marginal de \(\beta\) muestra con claridad que este parámetro no resultó significativo ya que el 0 está dentro de su intervalo de credibilidad
Inferencia sobre las posiciones latentes
library(MCMCpack)
# Transformación de Procrustes
B <- dim(modelo_ld$sample$Z)[1] # Número de muestras MCMC
n <- dim(modelo_ld$sample$Z)[2] # Número de vértices
d <- dim(modelo_ld$sample$Z)[3] # Dimensión latente
U0 <- scale(modelo_ld$mcmc.mle$Z, scale = TRUE, center = TRUE)
U.array <- array(data = NA, dim = c(B, n, d))
for (b in 1:B) {
U.array[b,,] <- MCMCpack::procrustes(
X = scale(modelo_ld$sample$Z[b,,], scale = TRUE, center = TRUE),
Xstar = U0,
translation = TRUE,
dilation = TRUE
)$X.new
}
U.pm <- apply(X = U.array, MARGIN = c(2, 3), FUN = mean)
# Colores (Asegurando que los valores estén en [0, 1])
rr <- atan2(U0[, 2], U0[, 1])
rr <- (rr - min(rr)) / (max(rr) - min(rr)) # Escalar a [0,1]
gg <- 1 - rr # Complemento
bb <- (U0[, 2]^2 + U0[, 1]^2)
bb <- (bb - min(bb)) / (max(bb)) # Escalar a [0,1]
aa <- 0.4 # Transparencia fija
# Adelgazamiento de la cadena
nthin <- 10
index_thin <- seq(from = nthin, to = B, by = nthin)
thinned_data <- do.call(rbind, lapply(index_thin, function(b) {
data.frame(
Dim1 = U.array[b,,1],
Dim2 = U.array[b,,2],
Vertex = factor(1:n),
Color = rgb(rr, gg, bb, alpha = aa)
)
}))
# Datos para las posiciones promedio
U_pm_df <- data.frame(
Dim1 = U.pm[, 1],
Dim2 = U.pm[, 2],
Vertex = 1:n,
Color = rgb(rr, gg, bb, alpha = 1)
)
# Primer panel: Con trayectorias y etiquetas
p1 <- ggplot() +
geom_point(data = thinned_data, aes(x = Dim1, y = Dim2, color = Color),
shape = 15, size = 0.8) + # Puntos de trayectorias
geom_text(data = U_pm_df, aes(x = Dim1, y = Dim2, label = Vertex),
color = "black", size = 3, fontface = "bold") + # Etiquetas en negro
scale_color_identity() +
labs(
x = "Dimensión 1",
y = "Dimensión 2",
title = "Posiciones latentes con trayectorias"
) +
theme_minimal()
# Segundo panel: Solo posiciones promedio con etiquetas pequeñas
p2 <- ggplot() +
geom_point(data = U_pm_df, aes(x = Dim1, y = Dim2, color = "blue"),
size = 3) + # Puntos tradicionales
geom_text(data = U_pm_df, aes(x = Dim1, y = Dim2, label = Vertex),
color = "black", size = 3, vjust = 1.5) + # Etiquetas pequeñas debajo de los puntos
scale_color_identity() +
labs(
x = "Dimensión 1",
y = "Dimensión 2",
title = "Posiciones promedio con etiquetas"
) +
theme_minimal()
# Combinar gráficos en dos paneles
grid.arrange(p1, p2, ncol = 2,
top = textGrob("Inferencia sobre las posiciones latentes",
gp = gpar(fontsize = 16, fontface = "bold") )
)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
Note que en la visualización de las posiciones latentes con trayectorias es difícil diferenciar el plano en grupos definidos, ya que las nubes de puntos se superponen entre si. Por su parte, en el gráfico de posiciones promedio si vemos a ciertos individuos muy cercanos de otros en el espacio latente.
# probabilidades de interacción (media posterior)
Pi <- matrix(0, n, n)
for (b in 1:B) {
bet <- modelo_ld$sample$beta[b]
for (i in 1:(n-1)) {
for (j in (i+1):n) {
lat <- sqrt(sum((modelo_ld$sample$Z[b,i,] - modelo_ld$sample$Z[b,j,])^2))
Pi[i,j] <- Pi[j,i] <- Pi[i,j] + expit(bet - lat)/B
}
}
}
# gráfico
rownames(Y) <- colnames(Y) <- 1:n
par(mfrow = c(1,2), oma = c(0, 0, 1, 0))
corrplot::corrplot(corr = Y,
type = "full",
col.lim = c(0,1),
method = "color",
tl.col = "black",
tl.pos = "n")
corrplot::corrplot(corr = Pi,
type = "full",
col.lim = c(0,1),
method = "color",
tl.col = "black",
tl.pos = "n")
mtext("Probabilidades de interacción (Real vs Modelo)",
side = 3, outer = TRUE, line = -1, cex = 1.3, font = 1)El modelo latente de distancia, a pesar de no contar con un parámetro \(\beta\) significativo, reconstruye bastante bien la matriz de adyacencia observada sin sobre ajustarse a esta. El modelo devuelve altas probabilidades de interacciones en la esquina superior izquierda de la matriz, la cual es la parte más conectada de la red de mafia.
# Hiperparámetros
a_sigma <- 2
b_sigma <- 1
a_tau <- 2
b_tau <- 1
# Ajustar el modelo usando Gibbs
n_iter <- 50000 + 10000
n_burn <- 10000
n_thin <- 5
#samples <- gibbs_sampler(Y, n_iter, n_burn, n_thin)
#save(samples, file = "Datos/samples_socio_mafia.RData")
load("Datos/samples_socio_mafia.RData")
round(mc_ee(samples$mu), 4)## var1
## 0.0047
## var1
## 0.0265
## var1
## 8e-04
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0071 0.0081 0.0092 0.0091 0.0101 0.0107
El modelo cuenta con errores estandar de Monte Carlo pequeños para los parámetros estimados.
# Calcular la log-verosimilitud para las muestras del muestreador
#log_lik <- log_likelihood(Y, samples)
#save(log_lik, file = "Datos/log_lik_mafia.RData")
load("Datos/log_lik_mafia.RData")
plot(log_lik, main = "Log-verosimilitud Red de consenso", xlab = "Iteración", ylab = "Log-verosimilitud", pch = 20, col = "cyan3")Por la estacionariedad de la cadena de máxima verosimilitud del modelo podemos afirmar que este convergió correctamente.
# Función para preparar los deltas
delta_dfs <- crear_delta_df(samples)
graficar_delta(delta_dfs, "Inferencia sobre los parámetros delta \n Red de consenso")## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
La red de mafia cuenta con varios parámetros de sociabilidad \(\delta_i\) que dentro de sus intervalos de credibilidad no incluyen el cero, por lo que para estos dan cuenta de una sociabilidad significativamente positiva.
Resalta el individuo 27, el cual cuyo IC de su parámetro se aleja más de los demás siendo un nodo más popular en la red de mafia. Este individuo 27 corresponde a ‘N18’ el cual hace parte de la cúpula ejecutiva de la familia Mistretta; y el que en análisis previos vimos que es el individuo con mayor grado en la red por una alta diferencia con los siguientes más conectados.
# Calcular la matriz promedio de theta_ij
#theta_avg <- compute_theta(samples)
#save(theta_avg, file = "Datos/theta_avg_mafia.RData")
load("Datos/theta_avg_mafia.RData")
graficar_theta_vs_y(theta_avg, Y, titulo = "Red de mafia")Las probabilidades de interacción resultantes del modelo de sociabilidad emulan la alta conectividad que se concentra en la parte inferior izquierda de la matriz de adyacencia original
# Calcular probabilidades de coagrupamiento
coclustering_probs <- compute_coclustering(samples)
# Estimar clusters usando Mclust
clusters <- estimate_clusters_mclust(coclustering_probs)
E <- compute_binary_E(clusters)
# Reordenar las matrices según los clusters estimados
coclustering_reordered <- reorder_matrix(coclustering_probs, clusters)
E_reordered <- reorder_matrix(E, clusters)
Y_reordered <- reorder_matrix(Y, clusters)
# Visualización de las tres matrices reordenadas
p1 <- plot_matrix(coclustering_reordered, "Co-clustering probabilities")
p2 <- plot_matrix(E_reordered, "Estimación puntual de clusters")
p3 <- plot_matrix(Y_reordered, "Matriz de adyacencia observada")
# Mostrar gráficos en un panel
gridExtra::grid.arrange(p1, p2, p3, ncol = 3)Al agrupar los \(\delta_i\) se obtienen las ya mencionadas probabilidades de co-agrupamiento. Se observa que las probabilidades posteriores de pertenencia al mismo cluster delimitan visualmente lo que en el gráfico dos son las estimación puntuales de los clusters. Resaltan los individuos que hacen parte del cluster 3 ya que dentro de las iteraciones del MCMC casi siempre resultaron en el mismo cluster juntos. La matriz de adyacencia ordenada a partir de esta partición muestra, por ejemplo, que los individuos del primer conglomerado están a duras penas conectados entre si.
# 1. Calcular las estimaciones y los intervalos de credibilidad
delta_mean <- colMeans(samples$delta) # Media posterior de delta
delta_ci95 <- apply(samples$delta, 2, quantile, probs = c(0.025, 0.975)) # IC al 95%
delta_df <- data.frame(Node = 1:n,
Delta_Est = delta_mean,
CI95_Lower = delta_ci95[1, ],
CI95_Upper = delta_ci95[2, ])
delta_df$Cluster <- as.factor(clusters) # Etiquetas de clusters como factores
delta_df <- delta_df[order(delta_df$Delta_Est), ]
delta_df$Order <- 1:n
ggplot(delta_df, aes(x = Order, y = Delta_Est, color = Cluster)) +
# IC
geom_segment(aes(x = Order, xend = Order,
y = CI95_Lower, yend = CI95_Upper),
size = 0.8) +
# Lineas pequeñas en los extremos
geom_segment(aes(x = Order - 0.2, xend = Order + 0.2,
y = CI95_Lower, yend = CI95_Lower),
size = 0.8) +
geom_segment(aes(x = Order - 0.2, xend = Order + 0.2,
y = CI95_Upper, yend = CI95_Upper),
size = 0.8) +
geom_point(size = 2) +
geom_text(aes(y = CI95_Upper + 0.1, label = Node),
size = 3,
hjust = 0.5) +
geom_hline(yintercept = 0,
linetype = "dashed",
color = "red") +
scale_color_manual(values = rainbow(length(unique(clusters)))) +
labs(title = "Inferencia sobre los parámetros delta",
x = NULL, # Eliminar etiquetas del eje x
y = expression(delta) ) +
theme_minimal(base_size = 14) +
theme(plot.title = element_text(hjust = 0.5, face = "bold"),
legend.position = "right",
axis.text.x = element_blank(),
axis.ticks.x = element_blank() )Los agrupaciones de los parámetros de sociabilidad generan 9 clases bien definidas de acuerdo al nivel de popularidad de los nodos. Tiene sentido que el cluster 1 se caracterice como la clase con los individuos menos conectados, comportamiento que vimos en la matriz de adyacencia ordenada.
# valores observados
dens_obs <- igraph::edge_density(graph = G, loops = F)
tran_obs <- igraph::transitivity(graph = G, type = "global")
asso_obs <- igraph::assortativity_nominal(graph = G,
types = as.integer(factor(vertex_attr(G, "familia"))),
directed = F)
mdis_obs <- igraph::mean_distance(graph = G, directed = F)
# Crear un dataframe con los valores observados
obs_values <- data.frame(
Estadistico = c("Densidad", "Transitividad", "Asortatividad", "DistanciaProm"),
Observado = c(dens_obs, tran_obs, asso_obs, mdis_obs)
)
# = = = = = = = = = = = = Función para la Bondad de ajuste = = = = = = = = = = = =
plot_gof <- function(stat_df, obs_values, title){
# Estadisticas en formato largo
stat_long <- stat_df %>%
pivot_longer(everything(),
names_to="Estadistico",
values_to="Valor")
# IC
ci_values <- stat_long %>%
group_by(Estadistico) %>%
summarise(CI_Lower = quantile(Valor,.025, na.rm = TRUE),
CI_Upper = quantile(Valor,.975, na.rm = TRUE) )
plot_df <- stat_long %>%
left_join(obs_values, by = "Estadistico") %>%
left_join(ci_values, by = "Estadistico")
# Grafico
ggplot(plot_df, aes(Valor))+
geom_histogram(aes(y = after_stat(density)),
bins = 20,
fill = "gray85",
color = "gray60")+
geom_vline(aes(xintercept = Observado),
colour = "red",
linewidth = 1)+
geom_vline(aes(xintercept = CI_Lower),
colour = "blue",
linetype = 2)+
geom_vline(aes(xintercept = CI_Upper),
colour = "blue",
linetype = 2)+
facet_wrap(~Estadistico, scales = "free",
labeller = labeller(Estadistico = c(DistanciaProm = "Distancia promedio",
Asortatividad = "Asortatividad: Familia/Clan")
) )+
labs(y = "Densidad",
title = title) +
theme_minimal() +
theme(plot.title = element_text( size = 16) )
}
tabla_gof <- function(stat_df, obs_values, caption){
# Resumen de las simulaciones
resumen <- tibble(Estadistico = names(stat_df),
Media = sapply(stat_df, mean, na.rm = TRUE),
SD = sapply(stat_df, sd, na.rm = TRUE),
IC_Inf = sapply(stat_df, quantile, probs = 0.025, na.rm = TRUE),
IC_Sup = sapply(stat_df, quantile, probs = 0.975, na.rm = TRUE))
# Unir con valores observados
tabla <- obs_values %>%
rename(Observado = Observado) %>%
left_join(resumen, by = "Estadistico")
tabla$Estadistico <- recode(tabla$Estadistico,
DistanciaProm = "Distancia promedio",
Asortatividad = "Asortatividad: Familia/Clan")
tabla <- tabla %>%
arrange(Estadistico)
tabla %>%
kbl(caption = caption, digits = 4, row.names = FALSE,
col.names = c("Estadistica", "Observado", "Media", "Desv. Est.", "IC Inferior", "IC Superior"),
align = c("l", rep("c", 5)) ) %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed"),
full_width = FALSE)
}
familia <- as.integer(factor(vertex_attr(G, "familia")))# simulación
B <- 5000
n <- vcount(G)
stat <- matrix(NA, B, 4)
set.seed(123)
for (b in 1:B) {
# grafo
g <- sample_gnp(n = n, p = theta_hat, directed = F, loops = F)
# estadísticos
stat[b,1] <- edge_density(g, loops = F)
stat[b,2] <- transitivity(g, type = "global")
stat[b,3] <- assortativity_nominal(g, types = familia, directed = F)
stat[b,4] <- mean_distance(g, directed = F)
}
# Dataframe con las estadísticas simuladas
stat_df <- data.frame(Densidad = stat[,1],
Transitividad = stat[,2],
Asortatividad = stat[,3],
DistanciaProm = stat[,4])
plot_gof(stat_df, obs_values,"Bondad de ajuste Erdos-Renyi")tabla_gof(stat_df = stat_df, obs_values, "Modelo 1: Resumen de las estadísticas observadas y simuladas con IC del 95%")| Estadistica | Observado | Media | Desv. Est. | IC Inferior | IC Superior |
|---|---|---|---|---|---|
| Asortatividad: Familia/Clan | 0.1573 | -0.0085 | 0.0330 | -0.0726 | 0.0565 |
| Densidad | 0.0320 | 0.0320 | 0.0017 | 0.0288 | 0.0355 |
| Distancia promedio | 4.1789 | 3.4109 | 0.1101 | 3.2151 | 3.6418 |
| Transitividad | 0.2785 | 0.0316 | 0.0079 | 0.0168 | 0.0478 |
El primer modelo que corresponde a un modelo de grafos aleatorios solo logra capturar la densidad de la red de mafia. Subestima la homofilia inducida por la familia/clan, el distanciamiento de los nodos y su capacidad para cerrar tríadas.
Como la probabilidad de interacción es constante e igual a la densidad, en promedio las redes simuladas tienen densidad igual a la observada (véase esto en la tabla).
# simulación
B <- 5000
n <- vcount(G)
stat <- matrix(NA, B, 4)
set.seed(123)
for (b in 1:B) {
# grafo
g <- sample_degseq(out.deg = grados, method = "vl")
# estadísticos
stat[b,1] <- edge_density(g, loops = F)
stat[b,2] <- transitivity(g, type = "global")
stat[b,3] <- assortativity_nominal(g, types = familia, directed = F)
stat[b,4] <- mean_distance(g, directed = F)
}
# Dataframe con las estadísticas simuladas
stat_df <- data.frame(Densidad = stat[,1],
Transitividad = stat[,2],
Asortatividad = stat[,3],
DistanciaProm = stat[,4])
plot_gof(stat_df, obs_values,"Bondad de ajuste Erdos-Renyi generalizado")tabla_gof(stat_df = stat_df, obs_values, "Modelo 2: Resumen de las estadísticas observadas y simuladas con IC del 95%")| Estadistica | Observado | Media | Desv. Est. | IC Inferior | IC Superior |
|---|---|---|---|---|---|
| Asortatividad: Familia/Clan | 0.1573 | -0.0342 | 0.0262 | -0.0828 | 0.0189 |
| Densidad | 0.0320 | 0.0320 | 0.0000 | 0.0320 | 0.0320 |
| Distancia promedio | 4.1789 | 3.1017 | 0.0456 | 3.0225 | 3.2024 |
| Transitividad | 0.2785 | 0.1200 | 0.0115 | 0.0973 | 0.1433 |
Dado que se fija una secuencia de grados, implícitamente se está fijando el número de aristas, por lo que todas las redes simuladas tienen la misma densidad que la red de consenso observada, aunque en el gráfico no se note este comportamiento.
Al igual que el modelo anterior, el modelo de grafos aleatorios generalizado solo captura de la densidad de la red por lo mencionado anteriormente. Subestima las otras características estructurales estudiadas por una diferencia considerable.
## Sampling ■■■■■■ 18% | ETA: 5sSampling ■■■■■■■ 21% | ETA: 5sSampling ■■■■■■■■
## 24% | ETA: 5sSampling ■■■■■■■■■ 26% | ETA: 5sSampling ■■■■■■■■■■ 29% | ETA:
## 5sSampling ■■■■■■■■■■ 32% | ETA: 5sSampling ■■■■■■■■■■■ 34% | ETA: 4sSampling
## ■■■■■■■■■■■■ 37% | ETA: 4sSampling ■■■■■■■■■■■■■ 40% | ETA: 4sSampling
## ■■■■■■■■■■■■■■ 43% | ETA: 4sSampling ■■■■■■■■■■■■■■■ 46% | ETA: 4sSampling
## ■■■■■■■■■■■■■■■ 48% | ETA: 4sSampling ■■■■■■■■■■■■■■■■ 51% | ETA: 3sSampling
## ■■■■■■■■■■■■■■■■■ 54% | ETA: 3sSampling ■■■■■■■■■■■■■■■■■■ 57% | ETA:
## 3sSampling ■■■■■■■■■■■■■■■■■■■ 59% | ETA: 3sSampling ■■■■■■■■■■■■■■■■■■■■ 62% |
## ETA: 3sSampling ■■■■■■■■■■■■■■■■■■■■ 65% | ETA: 3sSampling
## ■■■■■■■■■■■■■■■■■■■■■ 68% | ETA: 2sSampling ■■■■■■■■■■■■■■■■■■■■■■ 70% | ETA:
## 2sSampling ■■■■■■■■■■■■■■■■■■■■■■■ 73% | ETA: 2sSampling
## ■■■■■■■■■■■■■■■■■■■■■■■■ 76% | ETA: 2sSampling ■■■■■■■■■■■■■■■■■■■■■■■■■ 78% |
## ETA: 2sSampling ■■■■■■■■■■■■■■■■■■■■■■■■■ 81% | ETA: 1sSampling
## ■■■■■■■■■■■■■■■■■■■■■■■■■■ 84% | ETA: 1sSampling ■■■■■■■■■■■■■■■■■■■■■■■■■■■
## 87% | ETA: 1sSampling ■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 90% | ETA: 1sSampling
## ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 92% | ETA: 1sSampling
## ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 95% | ETA: 0sSampling
## ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 98% | ETA: 0s
stat <- matrix(NA, B, 4)
for(b in seq_len(B)){
g <- intergraph::asIgraph(simulated_networks[[b]])
# Calcular estadísticas
stat[b,1] <- edge_density(g, loops = F)
stat[b,2] <- transitivity(g, type = "global")
stat[b,3] <- assortativity_nominal(g, types = familia, directed = F)
stat[b,4] <- mean_distance(g, directed = F)
}
stat_df <- as.data.frame(stat)
colnames(stat_df) <- c("Densidad", "Transitividad", "Asortatividad", "DistanciaProm")
plot_gof(stat_df, obs_values,"Bondad de ajuste ERGM")tabla_gof(stat_df = stat_df, obs_values, "Modelo 3: Resumen de las estadísticas observadas y simuladas con IC del 95%")| Estadistica | Observado | Media | Desv. Est. | IC Inferior | IC Superior |
|---|---|---|---|---|---|
| Asortatividad: Familia/Clan | 0.1573 | 0.1260 | 0.0345 | 0.0584 | 0.1943 |
| Densidad | 0.0320 | 0.0320 | 0.0017 | 0.0288 | 0.0354 |
| Distancia promedio | 4.1789 | 3.4042 | 0.1066 | 3.2098 | 3.6271 |
| Transitividad | 0.2785 | 0.0760 | 0.0166 | 0.0460 | 0.1115 |
El modelo de grafos aleatorios exponencial además de la densidad logra capturar la homofilia inducida por la variable familia/clan a la que pertenecen los individuos de la red. Esto último era de esperarse ya que el modelo se construyó específicamente con los efectos homofílicos de la variable nodal mencionada. Note, por ejemplo, que la asortatividad media de los datos simulados (0.126) es cercana al valor observado (0.157).
No obstante este sigue subestimando la distancia geodésica y la transitividad de nuestra red.
Pi <- modelo_sbm$model_parameters[[Q]]$pi
Pi <- 0.5*(t(Pi) + Pi)
B <- 5000
n <- vcount(G)
stat <- matrix(NA, B, 4)
set.seed(123)
for (b in 1:B) {
# Red simulada
bs <- stats::rmultinom(n = 1, size = n, prob = alpha)
g_sim <- igraph::sample_sbm(n = n, pref.matrix = Pi, block.sizes = bs, directed = FALSE)
# estadísticos
stat[b,1] <- edge_density(g_sim, loops = F)
stat[b,2] <- transitivity(g_sim, type = "global")
stat[b,3] <- assortativity_nominal(g_sim, types = familia, directed = F)
stat[b,4] <- mean_distance(g_sim, directed = F)
}
# Dataframe con las estadísticas simuladas
stat_df <- data.frame(Densidad = stat[,1],
Transitividad = stat[,2],
Asortatividad = stat[,3],
DistanciaProm = stat[,4])
plot_gof(stat_df, obs_values,"Bondad de ajuste SBM")tabla_gof(stat_df = stat_df, obs_values, "Modelo 4: Resumen de las estadísticas observadas y simuladas con IC del 95%")| Estadistica | Observado | Media | Desv. Est. | IC Inferior | IC Superior |
|---|---|---|---|---|---|
| Asortatividad: Familia/Clan | 0.1573 | -0.0205 | 0.0306 | -0.0810 | 0.0410 |
| Densidad | 0.0320 | 0.0322 | 0.0060 | 0.0216 | 0.0450 |
| Distancia promedio | 4.1789 | 3.0501 | 0.2465 | 2.6879 | 3.6243 |
| Transitividad | 0.2785 | 0.1440 | 0.0226 | 0.0994 | 0.1883 |
Bao el criterio de bondad de ajuste, el modelo de bloques estocásticos se comporta igual que los modelos de la familia Erdos-Rényi debido a que solo logran capturar la densidad de la red observada. Las demás estadísticas estructurales las subestiman. Note así mismo que las distancias geodésicas promedio de las redes simuladas por el modelo son bastante dispersas por lo que su desviación estándar toma un valor alto de 0.24, y es lo que hace que se vea una cola a derecha.
# bondad de ajuste
B <- dim(modelo_ld$sample$Z)[1]
n <- dim(modelo_ld$sample$Z)[2]
d <- dim(modelo_ld$sample$Z)[3]
stat <- matrix(NA, B, 4)
set.seed(123)
for (b in 1:B) {
# intercepto
bet <- modelo_ld$sample$beta[b]
# simular datos
Ar <- matrix(0, n, n)
for (i in 1:(n-1)) {
for (j in (i+1):n){
lat <- sqrt(sum((modelo_ld$sample$Z[b,i,] - modelo_ld$sample$Z[b,j,])^2))
Ar[i,j] <- Ar[j,i] <- rbinom(n = 1, size = 1, prob = expit(bet - lat))
}
}
gr <- igraph::graph_from_adjacency_matrix(adjmatrix = Ar, mode = "undirected")
# calcular estadísticos
stat[b,1] <- igraph::edge_density(graph = gr, loops = F)
stat[b,2] <- igraph::transitivity(graph = gr, type = "global")
stat[b,3] <- igraph::assortativity_nominal(gr, types = familia, directed = F)
stat[b,4] <- igraph::mean_distance(graph = gr, directed = F)
}
save(stat, file = "Datos/stats_ldm_mafia.RData")load("Datos/stats_ldm_mafia.RData")
# Crear un dataframe con las estadísticas simuladas
stat_df <- data.frame(
Densidad = stat[,1],
Transitividad = stat[,2],
Asortatividad = stat[,3],
DistanciaProm = stat[,4])
plot_gof(stat_df, obs_values,"Bondad de ajuste LDM")tabla_gof(stat_df = stat_df, obs_values, "Modelo 5: Resumen de las estadísticas observadas y simuladas con IC del 95%")| Estadistica | Observado | Media | Desv. Est. | IC Inferior | IC Superior |
|---|---|---|---|---|---|
| Asortatividad: Familia/Clan | 0.1573 | 0.1032 | 0.0377 | 0.0319 | 0.1782 |
| Densidad | 0.0320 | 0.0320 | 0.0023 | 0.0277 | 0.0366 |
| Distancia promedio | 4.1789 | 3.7752 | 0.1844 | 3.4514 | 4.1673 |
| Transitividad | 0.2785 | 0.1333 | 0.0206 | 0.0938 | 0.1762 |
El modelo latente de distancia, que fue el mejor comportado para la red de consenso del punto anterior, en este caso resalta nuevamente por su desempeño sin cumplir todos los criterios. Junto con el ERGM, estos son los dos únicos modelos que tratan de emular la leve homofilia de la familia/clan de los individuos. Igualmente, los intervalos de credibilidad de la distancia geodésica promedio por un pelo y contienen al valor observado (limite superior de 4.16 vs. observado de 4.17). No obstante, se repite el patrón de subestimar la transitividad de la red.
B <- length(samples$mu)
n <- length(samples$delta[1, ])
stat <- matrix(NA, B, 4)
set.seed(123)
for(b in 1:B){
mu <- samples$mu[b]
delta <- samples$delta[b, ]
# Matriz de probabilidades
P <- pnorm(mu + outer(delta, delta, "+"))
# Simular red
Y_sim <- matrix(rbinom(length(P), 1, P), n, n)
Y_sim[lower.tri(Y_sim)] <- t(Y_sim)[lower.tri(Y_sim)]
diag(Y_sim) <- 0
g_sim <- graph_from_adjacency_matrix(Y_sim, mode = "undirected")
# Estadísticos
stat[b,1] <- edge_density(g_sim, loops = F)
stat[b,2] <- transitivity(g_sim, type = "global")
stat[b,3] <- assortativity_degree(g_sim, directed = F)
stat[b,4] <- mean_distance(g_sim, directed = F)
}
save(stat, file = "Datos/stat_socio_mafia.rData")load("Datos/stat_socio_mafia.rData")
stat_df <- as.data.frame(stat)
names(stat_df) <- c("Densidad", "Transitividad", "Asortatividad", "DistanciaProm")
plot_gof(stat_df, obs_values,"Bondad de ajuste Sociabilidad")tabla_gof(stat_df = stat_df, obs_values, "Modelo 6: Resumen de las estadísticas observadas y simuladas con IC del 95%")| Estadistica | Observado | Media | Desv. Est. | IC Inferior | IC Superior |
|---|---|---|---|---|---|
| Asortatividad: Familia/Clan | 0.1573 | -0.1956 | 0.0375 | -0.2681 | -0.1191 |
| Densidad | 0.0320 | 0.0321 | 0.0023 | 0.0278 | 0.0367 |
| Distancia promedio | 4.1789 | 2.8458 | 0.0925 | 2.6769 | 3.0388 |
| Transitividad | 0.2785 | 0.1096 | 0.0154 | 0.0806 | 0.1409 |
El último modelo, el modelo de sociabilidad, cierra con una bondad de ajuste deficiente. Este modelo solo logra capturar la densidad de la red de mafia de forma muy precisa, mientras que subestima todas las demás características.
| Modelo | Asort. Clan | Densidad | Dist. Promedio | Transitividad | |
|---|---|---|---|---|---|
| Erdos-Rényi | ❌ | ✅ | ❌ | ❌ | |
| Erdos-Rényi generalizado | ❌ | ✅ | ❌ | ❌ | |
| ERGM | ✅ | ✅ | ❌ | ❌ | |
| SBM | ❌ | ✅ | ❌ | ❌ | |
| LDM | ✅ | ✅ | - | ❌ | |
| Modelo de sociabildiad | ❌ | ✅ | ❌ | ❌ |
En general, ninguno de los modelos ajustados recuperó la transitividad de la red de mafia y solo el modelo latente de distancia (LDM) trató de capturar la estructura de la distancia. Dentro de los considerados, el mejor modelo en términos de bondad de ajuste vuelve a ser el LDM. Aunque es importante notar que los modelos deben mejorarse para caracterizar toda la estructura de la red.
Para evaluar la capacidad predictiva de los modelos ajustados se lleva a cabo el procedimiento de validación cruzada. Esta se correrá sobre los mismos \(L=5\) folds para hacer los análisis comparables.
## [1] 10153
# Folds
L <- 5
# conformación de folds
set.seed(123)
fold_index_vec <- sample(x = 1:L, size = M, replace = T)
fold_index_mat <- matrix(data = 0, nrow = n, ncol = n)
fold_index_mat[lower.tri(fold_index_mat)] <- fold_index_vec
# distribución de folds
(tab <- table(fold_index_vec))## fold_index_vec
## 1 2 3 4 5
## 2095 2084 2024 1921 2029
La separación de las 101053 posibles aristas totales en los 5 folds es bastante homogénea en cuanto a los tamaños de estos para 4 de los 5 folds, ya que en uno de estos se tienen en cuenta aproximadamente 1000 aristas menos que en los demás.
# viz de folds a través de la matriz de adyacencia
corrplot::corrplot(corr = fold_index_mat/L,
col.lim = c(0,1),
method = "color",
tl.pos = "n",
cl.pos = "n") ## 0 1
## [1,] 0.970 0.030
## [2,] 0.968 0.032
## [3,] 0.970 0.030
## [4,] 0.964 0.036
## [5,] 0.967 0.033
La repartición de los folds se llevo a cabo correctamente ya que la matriz no muestra agrupaciones de ningún tipo.
A continuación se obtienen las matrices de probabilidad de
interacción para cada fold para cada modelo. Estas se guardan
al final en una lista llamada modelos que se usará
posteriormente para graficar las curvas ROC y calcular el área bajo la
curva AUC.
IP1 <- vector(mode = "list", L)
IP2 <- vector(mode = "list", L)
B <- 500
set.seed(123)
for (l in 1:L) {
cat("Ajustando fold", l, "de", L, "...\n")
# Datos de entrenamiento
y_train <- y
y_train[fold_index_vec == l] <- NA
Y_train <- matrix(data = 0, nrow = n, ncol = n)
Y_train[lower.tri(Y_train)] <- y_train
Y_train <- Y_train + t(Y_train)
# ========================= Modelo Erdos-Renyi =========================
# ajuste del modelo
theta_hat <- mean(y_train, na.rm = T)
# predecir
n_test <- tab[l]
inter_prob <- rep(0, n_test)
for (b in 1:B)
inter_prob <- inter_prob + rbinom(n = n_test, size = 1, prob = theta_hat)/B
IP1[[l]] <- inter_prob
# ========================= Modelo Erdos-Renyi Generealziado =========================
# ajuste del modelo
sec <- rowSums(x = Y_train, na.rm = T)
# correcciones
sec[sec == 0] <- 1
if (sum(sec) %% 2 != 0) sec[1] <- sec[1] + 1
# predecir
n_test <- tab[l]
inter_prob <- rep(0, n_test)
for (b in 1:B) {
YY <- as.matrix(as_adjacency_matrix(sample_degseq(out.deg = sec, method = "vl")))
yy <- YY[lower.tri(YY)]
inter_prob <- inter_prob + yy[fold_index_vec == l]/B
}
IP2[[l]] <- inter_prob
}# ========================= Modelo ERGM =========================
B <- 500
library(parallel)
library(doParallel)
cl <- makeCluster(detectCores() - 1)
registerDoParallel(cl)
IP_ergm <- foreach(l = 1:L, .packages = c("network", "ergm")) %dopar% {
y_train <- y
y_train[fold_index_vec == l] <- NA
Y_train <- matrix(NA, n, n)
Y_train[lower.tri(Y_train)] <- y_train
Y_train[upper.tri(Y_train)] <- t(Y_train)[upper.tri(Y_train)]
diag(Y_train) <- 0
net_train <- network::as.network(Y_train, directed = FALSE, matrix.type = "adjacency")
network::set.vertex.attribute(net_train, attrname = "familia", value = familia)
# Intento 1: rápido, con init "caliente"
ergm_fit <- tryCatch({
ergm(net_train ~ edges + nodematch("familia", diff = TRUE),
control = control.ergm(
MCMC.samplesize = 2000,
MCMC.burnin = 5000,
MCMLE.maxit = 10
))
}, error = function(e) NULL)
# Intento 2 (fallback): más conservador, sin init fijo, más burn-in
if (is.null(ergm_fit)) {
ergm_fit <- tryCatch({
ergm(net_train ~ edges + nodematch("familia", diff = TRUE),
control = control.ergm(
MCMC.samplesize = 4000,
MCMC.burnin = 20000,
MCMC.interval = 2000,
MCMLE.maxit = 20
))
}, error = function(e) {
message("Fold ", l, " falló en ambos intentos: ", e$message)
NULL
})
}
n_test <- sum(fold_index_vec == l)
if (is.null(ergm_fit)) {
return(rep(NA, n_test))
}
sim_nets <- simulate(ergm_fit, nsim = B, output = "network")
inter_prob <- rep(0, n_test)
for (b in 1:B) {
YY <- as.matrix.network(sim_nets[[b]])
yy <- YY[lower.tri(YY)]
inter_prob <- inter_prob + yy[fold_index_vec == l] / B
}
inter_prob
}
stopCluster(cl)# ========================= Modelo de Bloques Estocásticos =========================
IP_sbm <- vector(mode = "list", L)
set.seed(123)
for (l in 1:L) {
# matriz de entrenamiento: díadas de prueba codificadas como AUSENCIA (0)
Y_train <- Y
Y_train[fold_index_mat == l] <- 0
Y_train[t(fold_index_mat) == l] <- 0
diag(Y_train) <- 0
modelo_sbm <- blockmodels::BM_bernoulli(membership_type = "SBM_sym",
adj = Y_train,
verbosity = 0,
plotting = "")
modelo_sbm$estimate()
Q_star <- which.max(modelo_sbm$ICL)
Z <- modelo_sbm$memberships[[Q_star]]$Z
Pi <- modelo_sbm$model_parameters[[Q_star]]$pi
Prob <- Z %*% Pi %*% t(Z)
prob_vec <- Prob[lower.tri(Prob)]
IP_sbm[[l]] <- prob_vec[fold_index_vec == l]
cat("SBM: fold", l, "completado — Q* =", Q_star, "\n")
}# ========================= Modelo Latente de Distancia =========================
IP_ld <- vector(mode = "list", L)
set.seed(123)
for (l in 1:L) {
y_train <- y
y_train[fold_index_vec == l] <- NA
Y_train <- matrix(NA, n, n)
Y_train[lower.tri(Y_train)] <- y_train
Y_train[upper.tri(Y_train)] <- t(Y_train)[upper.tri(Y_train)]
diag(Y_train) <- 0
net_train <- network::as.network(Y_train, directed = FALSE, matrix.type = "adjacency")
modelo_ld <- ergmm(net_train ~ euclidean(d = 2, G = 0),
seed = 123,
control = ergmm.control(burnin = 1000,
interval = 10,
sample.size = 5000))
Prob <- predict(modelo_ld) # matriz n x n de probabilidades estimadas
prob_vec <- Prob[lower.tri(Prob)]
IP_ld[[l]] <- prob_vec[fold_index_vec == l]
cat("M. Latente Distancia: fold", l, "completado\n")
}# Código modificado por eficiencia computacional
# ========================= Modelo de Sociabilidad =========================
gibbs_sampler_missing <- function(Y, n_iter, n_burn, n_thin, valid_indices, indices,
log_file = NULL, fold_id = NA) {
n <- nrow(Y)
mu <- 0
delta <- rnorm(n, 0, 1)
sigma2 <- 1
tau2 <- 1
z <- matrix(0, n, n)
# Pre-extraer índices UNA sola vez (evita indexar indices[idx,] cada iteración)
i_valid <- indices[valid_indices, 1]
j_valid <- indices[valid_indices, 2]
n_samples <- (n_iter - n_burn) / n_thin
prob_sum <- numeric(length(valid_indices)) # acumulador en vez de matriz completa
n_acc <- 0
log_fn <- function(msg) {
if (!is.null(log_file)) {
cat(msg, file = log_file, append = TRUE)
}
}
log_fn(sprintf("[Fold %s] Iniciando muestreador de Gibbs...\n", fold_id))
# Reporta progreso cada 5% en vez de cada 10%, y con timestamp
report_every <- max(1, round(n_iter / 20))
for (t in 1:n_iter) {
# --- Imputación vectorizada de aristas faltantes (no es sample_z) ---
eta <- mu + delta[i_valid] + delta[j_valid]
draws <- rbinom(length(i_valid), 1, pnorm(eta))
Y[cbind(i_valid, j_valid)] <- draws
Y[cbind(j_valid, i_valid)] <- draws
# --- Pasos del Gibbs (sin modificar) ---
z <- sample_z(Y, mu, delta, z)
mu <- sample_mu(z, delta, sigma2)
delta <- sample_delta(z, mu, tau2, delta)
sigma2 <- sample_sigma2(mu)
tau2 <- sample_tau2(delta)
if (t > n_burn && (t - n_burn) %% n_thin == 0) {
prob_sum <- prob_sum + Y[cbind(i_valid, j_valid)]
n_acc <- n_acc + 1
}
if (t %% report_every == 0) {
log_fn(sprintf("[Fold %s] %s | Progreso: %d%% (iter %d/%d)\n",
fold_id, format(Sys.time(), "%H:%M:%S"),
round(100 * t / n_iter), t, n_iter))
}
}
log_fn(sprintf("[Fold %s] Muestreador completado.\n", fold_id))
return(prob_sum / n_acc) # ya es la probabilidad media, no la matriz de muestras
}
# Hiperparámetros
a_sigma <- 2; b_sigma <- 1
a_tau <- 2; b_tau <- 1
set.seed(123)
# Folds
fold_index_mat_sym <- fold_index_mat + t(fold_index_mat)
indices <- which(upper.tri(Y, diag = FALSE), arr.ind = TRUE)
n_edges <- nrow(indices)
fold_assignment <- fold_index_mat_sym[indices]
folds <- lapply(1:L, function(l) which(fold_assignment == l))
# Carpeta de logs y resultados parciales
dir.create("Datos/logs", showWarnings = FALSE)
n_cores <- min(L, parallel::detectCores() - 1)
cl <- makeCluster(n_cores, outfile = "") # habilita stderr/stdout del cluster
registerDoParallel(cl)
clusterExport(cl,
varlist = c("Y", "indices", "folds", "gibbs_sampler_missing",
"sample_z", "sample_mu", "sample_delta",
"sample_sigma2", "sample_tau2",
"a_sigma", "b_sigma", "a_tau", "b_tau"),
envir = environment())
clusterEvalQ(cl, { library(pROC) })
results <- foreach(l = 1:L, .packages = "pROC") %dopar% {
log_file <- sprintf("Datos/logs/fold_%d_log.txt", l)
cat(sprintf("=== Fold %d iniciado: %s ===\n", l, Sys.time()),
file = log_file, append = FALSE)
valid_indices <- folds[[l]]
Y_train <- Y
for (idx in valid_indices) {
i <- indices[idx, 1]; j <- indices[idx, 2]
Y_train[i, j] <- NA
Y_train[j, i] <- NA
}
# Parámetros reducidos: ganancia de tiempo directa (ver estimación previa)
n_iter <- 20000
n_burn <- 5000
n_thin <- 5
prob_missing <- gibbs_sampler_missing(Y_train, n_iter, n_burn, n_thin,
valid_indices, indices,
log_file = log_file, fold_id = l)
Y_valid_values <- Y[indices[valid_indices, ]]
roc_curve <- roc(Y_valid_values, prob_missing)
roc_data <- data.frame(FPR = 1 - roc_curve$specificities,
TPR = roc_curve$sensitivities)
fold_result <- list(auc = as.numeric(auc(roc_curve)), roc = roc_data)
# Checkpoint: se guarda apenas termina, sin esperar a los demás folds
saveRDS(fold_result, file = sprintf("Datos/fold_%d_result.rds", l))
cat(sprintf("=== Fold %d completado: %s ===\n", l, Sys.time()),
file = log_file, append = TRUE)
fold_result
}
stopCluster(cl)
auc_values <- sapply(results, `[[`, "auc")
roc_list <- lapply(results, `[[`, "roc")
save(auc_values, file = "Datos/auc_values_mafia.RData")
save(roc_list, file = "Datos/roc_list_mafia.RData")Obtenidos los insumos necesarios, a continuación se obtienen las curvas ROC para los modelos ajustados. Esta curva es una representación gráfica de la relación entre la Tasa de verdadderos positivos y la Tasa de falsos positivos.
modelos <- list("Erdos-Renyi" = IP1,
"Erdos-Renyi Generalizado" = IP2,
"ERGM" = IP_ergm,
"SBM" = IP_sbm,
"LDM" = IP_ld)
save(modelos, file = "Datos/modelos_val_cru_mafia.RData")load("Datos/modelos_val_cru_mafia.RData")
plots <- list()
for (nombre in names(modelos)) {
IP <- modelos[[nombre]]
roc_df <- data.frame()
aucs <- c()
for (l in 1:L) {
if (any(is.na(IP[[l]]))) {message(nombre, ": fold ", l, " omitido")
next}
y_test <- y[fold_index_vec == l]
pred <- ROCR::prediction(IP[[l]], y_test)
perf <- ROCR::performance(pred, "tpr", "fpr")
roc_df <- rbind(roc_df,
data.frame(FPR = perf@x.values[[1]],
TPR = perf@y.values[[1]],
Fold = factor(l) ) )
aucs <- c(aucs, ROCR::performance(pred, "auc")@y.values[[1]])
}
auc_mean <- mean(aucs)
n_folds <- length(aucs)
plots[[nombre]] <-
ggplot(roc_df,
aes(FPR, TPR, colour = Fold)) +
geom_line(linewidth = 1, alpha = .8) +
geom_abline(slope = 1,
intercept = 0,
linetype = "dashed",
colour = "grey50") +
coord_equal() +
annotate("label", x = .8, y = .12,
label = paste0("AUC = ", round(auc_mean,4), "\n", n_folds,"/",L," folds"),
size = 4) +
labs(title = nombre,
x = "Tasa de falsos positivos",
y = "Tasa de verdaderos positivos",
colour = "Fold") +
theme_minimal(base_size = 14) +
theme(plot.title = element_text(face = "bold", hjust = .5),
legend.position = "bottom",
panel.grid.minor = element_blank() )
}roc_df <- do.call(rbind,
lapply(1:L, function(l)
cbind(roc_list[[l]], Fold = factor(l))) )
auc_mean <- mean(auc_values)
n_folds <- length(auc_values)
plots[["Modelo de Sociabilidad"]] <-
ggplot(roc_df,
aes(x = FPR, y = TPR, colour = Fold)) +
geom_line(linewidth = 1, alpha = 0.8) +
geom_abline(slope = 1,
intercept = 0,
linetype = "dashed",
colour = "grey50") +
coord_equal() +
annotate("label", x = 0.72, y = 0.12,
label = paste0("AUC = ", round(auc_mean, 4), "\n", n_folds, "/", L, " folds"),
size = 4) +
labs(title = "Modelo de Sociabilidad",
x = "Tasa de falsos positivos",
y = "Tasa de verdaderos positivos",
colour = "Fold") +
theme_minimal(base_size = 14) +
theme(plot.title = element_text(face = "bold", hjust = 0.5),
legend.position = "bottom",
panel.grid.minor = element_blank() )
grid.arrange(grobs = plots, ncol = 3)## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database
## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database
## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database
## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database
## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database
## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database
A partir de las curvas de ROC obtendias se observa una clara diferencia en la capacidad predictiva de los modelos.
Empecemos el análisis con los modelos Erdos-Rényi y ERGM, los cuales se estarían quedando con el reconocimiento a los modelos más precarios en términos predictivos. Sus AUC promedio (área bajo la curva) toman valores de 0.48 y 0.6 respectivamente, los cuales al ser tan cercanos a 0.5 ponen en evidencia su nula capacidad discriminatoria equivalente a clasificar casos positivos y negativos al azar. En otras palabras, los modelos no logran establecer un buen compromiso entre detectar enlaces reales y evitar predecir enlaces inexistentes.
Dado que la tasa de falsos negativos (TFN) = 1− tasa de verdaderos positivos (TPR), una curva que permanece cerca de la esquina superior izquierda indica simultáneamente una baja tasa de falsos positivos y una baja tasa de falsos negativos.
El LDM muestra una mejora importante. En los primeros incrementos de la tasa de falsos positivos ya consigue sensibilidades relativamente altas. Sin embargo, para alcanzar sensibilidades cercanas al 100 % la tasa de falsos positivos aumenta de manera apreciable.
Los demás modelos (modelo de sociabilidad, modelo Erdos-Rényi generalizado y modelo de bloques estocásticos) son los mejores en capacidad predictiva, con comportamientos consistentes entre los cinco folds lo que sugiere que su desempeño es poco sensible a la partición utilizada en la validación cruzada.
En conjunto, los resultados muestran que los modelos que incorporan heterogeneidad entre nodos o estructura comunitaria logran una mejor discriminación entre enlaces y no enlaces.