Script para Social Network Analysis - Live Coding no canal EstaTiDados - 2020-09-16

estatidados_sna.R

############################
#     Live EstaTiDados     #
#     Neylson Crepalde     #
#         REDES !!!!       #
############################


library(tidyverse)
library(sna)
library(igraph)
library(ggplot2)
library(blockmodeling)
devtools::install_github("aslez/concoR")
library(concoR)

#### Vamos brincar com os dados de advogados de Lazega
#### https://www.stats.ox.ac.uk/~snijders/siena/Lazega_lawyers_data.htm

# Fazendo download e unzip dos arquivos
url = "https://www.stats.ox.ac.uk/~snijders/siena/LazegaLawyers.zip"
temp = tempfile()
download.file(url, temp)
unzip(temp, exdir = "./data")

# Lê a matriz de adjacências 
le_rede = function(file) {
  admat = read_table(file, col_names = F)
  admat = as.matrix(admat)
  rownames(admat) = colnames(admat)
  return(admat)
}

work = le_rede("data/ELwork.dat")
adv = le_rede("data/ELadv.dat")
frd = le_rede("data/ELfriend.dat")

# Transforma em objetos IGRAPH
gwork = graph_from_adjacency_matrix(work)
gadv  = graph_from_adjacency_matrix(adv)
gfrd  = graph_from_adjacency_matrix(frd)

# Visualizando
myplot = function(g, title="") {
  plot(g, vertex.size=6, vertex.label = NA,
       vertex.color = adjustcolor('red', .6),
       edge.color = adjustcolor('grey', .8),
       edge.arrow.size=.3,
       main=title)
}

myplot(gwork, "Work net")
myplot(gadv, "Advice net")
myplot(gfrd, "Friendship net")


# Lê os atributos
atributos = read_table("data/ELattr.dat", 
                       col_names = c(
                         "seniority", "status", "gender", "office",
                         "years", "age", "practice", "lschool"
                       ),
                       col_types = cols(
                         seniority = col_integer(),
                         status    = col_character(),
                         gender    = col_character(),
                         office    = col_character(),
                         years     = col_integer(),
                         age       = col_integer(),
                         practice  = col_character(),
                         lschool   = col_character()
                       ))

# Vincula os atributos a cada rede
vincula_atributos = function(g) {
  V(g)$seniority = atributos$seniority
  V(g)$status    = atributos$status
  V(g)$gender    = atributos$gender
  V(g)$office    = atributos$office
  V(g)$years     = atributos$years
  V(g)$age       = atributos$age
  V(g)$practice  = atributos$practice
  V(g)$lschool   = atributos$lschool
  return(g)
}

lista_novos = lapply(list(gwork, gadv, gfrd), vincula_atributos)
gwork = lista_novos[[1]]
gadv  = lista_novos[[2]]
gfrd  = lista_novos[[3]]

###################################################
# Investigando a rede WORK
gwork
diameter(gwork)
mean(degree(gwork))
edge_density(gwork)

ggplot() +
  geom_histogram(aes(degree(gwork)), color='white', bins=15) +
  labs(title='Distribuição de grau - Work', x='Grau')

sort(degree(gwork), dec=T)

ggplot() +
  geom_histogram(aes(betweenness(gwork)), color='white', bins=15) +
  labs(title='Distribuição de betweeness - Work', x='Grau')

sort(betweenness(gwork), dec=T)

ggplot() +
  geom_histogram(aes(constraint(gwork)), color='white', bins=15) +
  labs(title='Distribuição de Constraint - Work', x='Grau')

sort(constraint(gwork), dec=F)


## Análise de cluster
ceb = cluster_edge_betweenness(gwork)
cle = cluster_leading_eigen(gwork)
cwt = cluster_walktrap(gwork)

?plot.igraph
plot(gwork, vertex.size=6, vertex.label = NA,
     vertex.color = adjustcolor('red', .6),
     edge.color = adjustcolor('grey', .8),
     edge.arrow.size=.3,
     mark.groups = cwt,
     main="Work Network")

############################################################
# Identificação de Equivalência Estrutural com BLOCKMODELING
mat = as.matrix(get.adjacency(gwork))


class4 <- optRandomParC(M=mat, k=4, rep=10, approach="ss", blocks="com")
class5 <- optRandomParC(M=mat, k=5, rep=10, approach="ss", blocks="com")
class6 <- optRandomParC(M=mat, k=6, rep=10, approach="ss", blocks="com")


par(mfrow=c(1,3)) # set the plot window for one row and two columns
plot(class4, main="")
title("4 Block Partition")
plot(class5, main="") 
title("5 Block Partition")
plot(class6, main="") 
title("6 Block Partition")
par(mfrow=c(1,1)) # reset the plot window back to one row and one column

l = layout_with_fr(gwork)

par(mfrow=c(1,3)) # set the plot window for one row and two columns
plot(gwork, vertex.size=6, vertex.label = NA,
     vertex.color = class4$best$best1$clu,
     edge.color = adjustcolor('grey', .8),
     edge.arrow.size=.3,
     #mark.groups = cwt,
     layout = l,
     main="Blockmodeling - 4 groups")
plot(gwork, vertex.size=6, vertex.label = NA,
     vertex.color = class5$best$best1$clu,
     edge.color = adjustcolor('grey', .8),
     edge.arrow.size=.3,
     #mark.groups = cwt,
     layout = l,
     main="Blockmodeling - 5 groups")
plot(gwork, vertex.size=6, vertex.label = NA,
     vertex.color = class6$best$best1$clu,
     edge.color = adjustcolor('grey', .8),
     edge.arrow.size=.3,
     #mark.groups = cwt,
     layout = l,
     main="Blockmodeling - 6 groups")
par(mfrow=c(1,1))


# CONCOR - Convergence of Iterated Correlations
gwork_pc = delete_vertices(gwork, V(gwork)[degree(gwork) == 0]) # deleta um isolado
m0 = as.matrix(get.adjacency(gwork_pc))

bm = concor_hca(list(m0), max.iter = 50, p=2)
partitions = sna::blockmodel(m0, ec = bm$block)
plot(partitions)


#### EXTRA - Prova de Convergência CONCOR
gwork_pc = delete_vertices(gwork, V(gwork)[degree(gwork) == 0])
m0 = as.matrix(get.adjacency(gwork_pc))

m_new = matrix(nrow=nrow(m0), ncol=ncol(m0))
for (i in 1:15){
  if (i == 1) {
    for (i in 1:ncol(m0)) {
      for (j in 1:ncol(m0)) {
        m_new[i,j] = cor(m0[,i], m0[,j])
      }
    }
  } else {
    for (i in 1:ncol(m0)) {
      for (j in 1:ncol(m0)) {
        m_new[i,j] = cor(m_new[,i], m_new[,j])
      }
    }
  }
}

print(m_new)
####################################################################
# Testes de Hipótese 
# ERGMS - Exponential Random Graph Models
# Cenas do próximo episódio...