Librerias requeridas

igraph, ggplot2, knitr, stringr, dplyr, readr, devtools

Funciones

Datos

Los datos originales fueron obtenidos de enlace. Posteriormente, fueron depurados y organizados en un archivo .txt disponible en el siguiente enlace.

url <-'https://www.farmaciayestadistica.com/public/files/bases/resistance.txt'
data <- read.csv(url, header=TRUE, sep = "\t", 
                 colClasses = c("character", "numeric", "numeric", "numeric", "numeric", 
                                "numeric", "numeric", "numeric", "numeric"))
tab <- data[,2:9]
rownames(tab) <- data[,1]

Visualización de la tabla de datos a travéz de un gráfico de calor

farben <- colorRampPalette(c("white", "blue"))
pal_blue <- farben(15)
heatmap(as.matrix(tab), 
        Colv = NA, Rowv = NA,    
        col = pal_blue,    
        scale = "none",          
        main = "Matriz de calor",
        cexRow = 0.8,
        cexCol = 0.8)

Generación de grafo bipartito a partir de la matriz de incidencias

bn <- graph_from_biadjacency_matrix(tab, weighted = TRUE)   
class(bn)  #Confirmación que es un grafo
## [1] "igraph"
V(bn)$type
##  [1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [13] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [25] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [37] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [49] FALSE FALSE FALSE FALSE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE

Visualización del grafo

shapes <- c("circle","square")
colors <- c("green","purple")
set.seed(123)
plot(bn,vertex.color=colors[V(bn)$type+1],
vertex.shape=shapes[V(bn)$type+1],
vertex.size=7,vertex.label.degree=-pi/2,
vertex.label.dist=0,vertex.label.cex=1,
edge.label = E(bn)$weight,
vertex.label.color="black",
layout=layout_with_kk(bn))

edge_density(bn)
## [1] 0.2248588
desctop <- function(df,ntv) {
mat <- matrix(0, nrow = ntv, ncol = ncol(df))
for (j in 1:ncol(df)) {
  sorted_indices <- order(df[, j], decreasing = TRUE)
  top_rows <- rownames(df)[sorted_indices[1:ntv]]
  mat[, j] <- top_rows
}
colnames(mat) <- colnames(df)
print(mat)
}
vec1 <- unique(as.vector(desctop(tab,3)))
##      X3GCRKP              FREC    CRPA              MRSA       
## [1,] "Russian_Federation" "India" "Romania"         "Indonesia"
## [2,] "India"              "Peru"  "Slovak_Republic" "China"    
## [3,] "Bulgaria"           "China" "Argentina"       "Romania"  
##      X3GCREC              PRSP      VRE             CRKP    
## [1,] "India"              "Turkey"  "United_States" "Greece"
## [2,] "Russian_Federation" "India"   "Ireland"       "India" 
## [3,] "Mexico"             "Romania" "Peru"          "Italy"
vec2 <- colnames(tab)
vec <- c(vec1,vec2)
ig <- induced_subgraph(graph = bn, vids = vec)
shapes <- c("circle","square")
colors <- c("green","purple")

set.seed(123)
plot(ig,vertex.color=colors[V(ig)$type+1],
vertex.shape=shapes[V(ig)$type+1],
vertex.size=7,vertex.label.degree=-pi/2,
vertex.label.dist=0,vertex.label.cex=1,
edge.label = E(ig)$weight,
vertex.label.color="black",
layout=layout_with_kk(ig))

edge_density(ig)
## [1] 0.4743083
ig <- induced_subgraph(graph = bn, vids = vec)
tam1 <- sqrt(rowSums(tab[vec1,]))
tam2 <- sqrt(colSums(tab[vec1,vec2]))
tam <- c(tam1,tam2)

shapes <- c("circle","square")
colors <- c(adjustcolor("green", alpha.f = 0.5),adjustcolor("purple", alpha.f = 0.5))
set.seed(123)
plot(ig,vertex.color=colors[V(ig)$type+1],
vertex.shape=shapes[V(ig)$type+1],
vertex.size=tam,vertex.label.degree=-pi/2,
vertex.label.dist=0,vertex.label.cex=1,
vertex.label.color="black",
layout=layout_with_kk(ig))

tam1 <- sqrt(rowSums(tab))
tam2 <- sqrt(colSums(tab))
tam <- c(tam1,tam2)

shapes <- c("circle","square")
colors <- c(adjustcolor("green", alpha.f = 0.5),adjustcolor("purple", alpha.f = 0.1))
set.seed(123)
plot(bn,vertex.color=colors[V(bn)$type+1],
vertex.shape=shapes[V(bn)$type+1],
vertex.size=tam,vertex.label.degree=-pi/2,
vertex.label.dist=0,vertex.label.cex=1,
vertex.label.color="black",
layout=layout_with_kk(bn))

Generación de los grafos individuales por país o por antibiótico

bn.pr <- bipartite_projection(bn)   
bn.country <- bn.pr$proj1   #Extracción del grafo de los paises
bn.antibiotic <- bn.pr$proj2    #Extracción del grafo de los antibioticos
#Mirar que las densidades don las mismas
edge_density(bn.pr$proj1)  #graph.density(bn.pr$proj1)   densidad tomada de la lista
## [1] 1
edge_density(bn.country)   #densidad tomada del grafo extraido
## [1] 1
edge_density(bn.pr$proj2)  #graph.density(bn.pr$proj1)   densidad tomada de la lista
## [1] 1
edge_density(bn.antibiotic)   #densidad tomada del grafo extraido
## [1] 1

Generación de las matrices de adyacencia individuales de los grafos individuales

Gráficas de cada uno de los grafos

shapes <- c("circle","square")
colors <- c("green","purple")

par(mfrow=c(1,2), mai = c(0.2, 0.2, 0.2, 0.2))

plot(bn.country,vertex.color="green",
vertex.shape="circle",main="Paises",
edge.label=E(bn.country)$weight,
vertex.size=15,vertex.label.degree=-pi/2,
vertex.label.dist=2,vertex.label.cex=0.8,
vertex.label.color="black")

plot(bn.antibiotic,vertex.color="purple",
vertex.shape="square",main="Antibioticos",
edge.label=E(bn.antibiotic)$weight,
vertex.size=15,vertex.label.degree=-pi/2,
vertex.label.dist=2,vertex.label.cex=0.8,
vertex.label.color="black")

Binarización de la matriz de biadyacencia

vu=40 #valor umbral
TA <- ifelse(tab >= vu, 1, 0)
Mady_country <- TA%*%t(TA)
Mady_antibiotic <- t(TA)%*%TA
g_country <- graph_from_adjacency_matrix(Mady_country, weighted = TRUE, 
                                            mode = "undirected", diag=TRUE) 
# Descripción del grafo
desc(g_country)
## El grafo no es dirigido, ponderado, , tiene  52 vertices y 405 aristas
is_simple(g_country)  #Ayuda a determinar si los grafos tienen bucles consigo mismos: aristas que conectan un vértice consigo mismo
## [1] FALSE
g_countrys <- simplify(g_country)
is_simple(g_countrys)
## [1] TRUE
###########
# Eliminar vértices aislados
g_countrys_sva <- delete.vertices(g_countrys, which(degree(g_countrys) == 0))  
## Warning: `delete.vertices()` was deprecated in igraph 2.0.0.
## ℹ Please use `delete_vertices()` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
# Plotear el grafo sin vértices aislados
plot(g_countrys_sva, vertex.label = V(g_countrys_sva)$name, 
     edge.label = E(g_countrys_sva)$weight)

# Clusterización
kc <- cluster_fast_greedy(g_countrys_sva)
str(kc)
## Class 'communities'  hidden list of 6
##  $ merges    : num [1:29, 1:2] 27 10 12 20 21 25 26 37 6 13 ...
##  $ modularity: num [1:30] -0.0383 -0.0361 -0.0337 -0.0301 -0.0252 ...
##  $ membership: num [1:30] 2 2 1 1 2 2 1 1 1 2 ...
##  $ names     : chr [1:30] "India" "China" "Russian_Federation" "Romania" ...
##  $ algorithm : chr "fast greedy"
##  $ vcount    : num 30
kc_fast_greedy <- cluster_fast_greedy(g_countrys_sva)
kc_leading_eigen <- cluster_leading_eigen(g_countrys_sva)
kc_walktrap <- cluster_walktrap(g_countrys_sva)
kc_louvain <- cluster_louvain(g_countrys_sva)
kc_label_prop <- cluster_label_prop(g_countrys_sva)
kc_spinglass <- cluster_spinglass(g_countrys_sva)
kc_optimal <- cluster_optimal(g_countrys_sva)
kc_infomap <- cluster_infomap(g_countrys_sva)

# gráficos
igraph_options(vertex.size = 10, vertex.frame.color = "black")
par(mfrow = c(3,3))
l = layout_with_dh(g_countrys_sva)
plot(g_countrys_sva, vertex.label = NA, layout = l, vertex.color = kc_fast_greedy$membership, main = paste0("fast greedy: ", "Mod = ", round(modularity(kc_fast_greedy), 4)))
plot(g_countrys_sva, vertex.label = NA, layout = l, vertex.color = kc_leading_eigen$membership, main = paste0("leading eigen: ", "Mod = ", round(modularity(kc_leading_eigen), 4)))
plot(g_countrys_sva, vertex.label = NA, layout = l, vertex.color = kc_walktrap$membership, main = paste0("walktrap: ", "Mod = ", round(modularity(kc_walktrap), 4)))
plot(g_countrys_sva, vertex.label = NA, layout = l, vertex.color = kc_louvain$membership, main = paste0("louvain: ", "Mod = ", round(modularity(kc_louvain), 4)))
plot(g_countrys_sva, vertex.label = NA, layout = l, vertex.color = kc_label_prop$membership, main = paste0("label prop: ", "Mod = ", round(modularity(kc_label_prop), 4)))
plot(g_countrys_sva, vertex.label = NA, layout = l, vertex.color = kc_spinglass$membership, main = paste0("spinglass: ", "Mod = ", round(modularity(kc_spinglass), 4)))
plot(g_countrys_sva, vertex.label = NA, layout = l, vertex.color = kc_optimal$membership, main = paste0("optimal: ", "Mod = ", round(modularity(kc_optimal), 4)))
plot(g_countrys_sva, vertex.label = NA, layout = l, vertex.color = kc_infomap$membership, main = paste0("infomap: ", "Mod = ", round(modularity(kc_infomap), 4)))

set.seed(1)
plot(x = kc, y = g_countrys_sva, vertex.size = 12)

set.seed(1)
plot(x = kc, y = g_countrys_sva, mark.groups = NULL, edge.color = "darkgray", vertex.size = 12)

# Generación de grafo a partir de matriz de adyacencia
g_country <- graph_from_adjacency_matrix(Mady_country, weighted = TRUE)
plot(g_country,vertex.color="purple",
vertex.shape="square",main="Antibioticos",
edge.label=E(g_country)$weight,
vertex.size=15,vertex.label.degree=-pi/2,
vertex.label.dist=2,vertex.label.cex=0.8,
vertex.label.color="black")

g_antibiotic <- graph_from_adjacency_matrix(Mady_antibiotic, weighted = TRUE, 
                                            mode = "undirected", diag=FALSE) 
# Descripción del grafo
desc(g_antibiotic)
## El grafo no es dirigido, ponderado, , tiene  8 vertices y 21 aristas
is_simple(g_antibiotic)  #Ayuda a determinar si los grafos tienen bucles consigo mismos: aristas que conectan un vértice consigo mismo
## [1] TRUE
g_antibiotics <- simplify(g_antibiotic)
is_simple(g_countrys)
## [1] TRUE
###########
# Eliminar vértices aislados
g_antibiotics_sva <- delete.vertices(g_antibiotics, which(degree(g_antibiotics) == 0))  

# Plotear el grafo sin vértices aislados
plot(g_antibiotics_sva, vertex.label = V(g_antibiotics_sva)$name, 
     edge.label = E(g_antibiotics_sva)$weight)

# Clusterización
kc <- cluster_fast_greedy(g_antibiotics_sva)
str(kc)
## Class 'communities'  hidden list of 6
##  $ merges    : num [1:6, 1:2] 3 4 8 10 6 12 1 2 9 5 ...
##  $ modularity: num [1:7] -0.17898 -0.12987 -0.09229 -0.04103 -0.00164 ...
##  $ membership: num [1:7] 2 2 2 2 2 1 1
##  $ names     : chr [1:7] "X3GCRKP" "FREC" "CRPA" "MRSA" ...
##  $ algorithm : chr "fast greedy"
##  $ vcount    : num 7
kc_fast_greedy <- cluster_fast_greedy(g_antibiotics_sva)
kc_leading_eigen <- cluster_leading_eigen(g_antibiotics_sva)
kc_walktrap <- cluster_walktrap(g_antibiotics_sva)
kc_louvain <- cluster_louvain(g_antibiotics_sva)
kc_label_prop <- cluster_label_prop(g_antibiotics_sva)
kc_spinglass <- cluster_spinglass(g_antibiotics_sva)
kc_optimal <- cluster_optimal(g_antibiotics_sva)
kc_infomap <- cluster_infomap(g_antibiotics_sva)

# gráficos
l = layout_with_dh(g_antibiotics_sva)
igraph_options(vertex.size = 30, vertex.frame.color = "black", vertex.label.cex = 1, vertex.label.color="black")
par(mfrow = c(3,3))
set.seed(123)
plot(g_antibiotics_sva, layout = l, vertex.color = kc_fast_greedy$membership, main = paste0("fast greedy: ", "Mod = ", round(modularity(kc_fast_greedy), 4)))
plot(g_antibiotics_sva, layout = l, vertex.color = kc_leading_eigen$membership, main = paste0("leading eigen: ", "Mod = ", round(modularity(kc_leading_eigen), 4)))
plot(g_antibiotics_sva, layout = l, vertex.color = kc_walktrap$membership, main = paste0("walktrap: ", "Mod = ", round(modularity(kc_walktrap), 4)))
plot(g_antibiotics_sva, layout = l, vertex.color = kc_louvain$membership, main = paste0("louvain: ", "Mod = ", round(modularity(kc_louvain), 4)))
plot(g_antibiotics_sva, layout = l, vertex.color = kc_label_prop$membership, main = paste0("label prop: ", "Mod = ", round(modularity(kc_label_prop), 4)))
plot(g_antibiotics_sva, layout = l, vertex.color = kc_spinglass$membership, main = paste0("spinglass: ", "Mod = ", round(modularity(kc_spinglass), 4)))
plot(g_antibiotics_sva, layout = l, vertex.color = kc_optimal$membership, main = paste0("optimal: ", "Mod = ", round(modularity(kc_optimal), 4)))
plot(g_antibiotics_sva, layout = l, vertex.color = kc_infomap$membership, main = paste0("infomap: ", "Mod = ", round(modularity(kc_infomap), 4)))

library(networkD3)
## Warning: package 'networkD3' was built under R version 4.3.3
t_antibiotics_sva <- get.data.frame(g_antibiotics_sva)
## Warning: `get.data.frame()` was deprecated in igraph 2.0.0.
## ℹ Please use `as_data_frame()` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
simpleNetwork(t_antibiotics_sva)
net_antibiotics <- simpleNetwork(t_antibiotics_sva)
saveNetwork(net_antibiotics,file = 'net_antibiotics.html', selfcontained=TRUE)
library(htmltools)
includeHTML("C:/Users/johnj/Documents/FARMACIA_Y_ESTADISTICA/05_Publicaciones/020_anagraf/020_anagraf_03/net_antibiotics.html")
## Warning: `includeHTML()` was provided a `path` that appears to be a complete HTML document.
## ✖ Path: C:/Users/johnj/Documents/FARMACIA_Y_ESTADISTICA/05_Publicaciones/020_anagraf/020_anagraf_03/net_antibiotics.html
## ℹ Use `tags$iframe()` to include an HTML document. You can either ensure `path` is accessible in your app or document (see e.g. `shiny::addResourcePath()`) and pass the relative path to the `src` argument. Or you can read the contents of `path` and pass the contents to `srcdoc`.
forceNetwork 
plot(g_antibiotic,vertex.color="purple",
vertex.shape="square",main="Antibioticos",
edge.label=E(g_antibiotic)$weight,
vertex.size=15,vertex.label.degree=-pi/2,
vertex.label.dist=2,vertex.label.cex=0.8,
vertex.label.color="black")

# Crear el grafo a partir de la matriz de adyacencia
#grafo <- graph.adjacency(Mady_antibiotic, weighted = TRUE, mode = "undirected", diag = FALSE)

# Plotear el grafo
#plot(grafo, vertex.label = colnames(matriz_adyacencia), vz=E(grafo)$weight)

Referencias

https://www.oecd-ilibrary.org/sites/9789264307599-6-en/index.html?itemId=/content/component/9789264307599-6-en

https://www.cdc.gov/drugresistance/across-the-world.html Kolaczyk, E. D., & Cs´ardi, G. (2014). Statistical analysis of network data with R (Vol. 65). New York: Springer.