#*******************************************************************************
# Analysis of Social Networks
# Introduction to network modeling
# V. Amati, A. Lomi

# This script illustrates how to generate networks from stylized and 
# simple network models and perform CUG test
#*******************************************************************************


# ------------------------------ WD and packages ----------------------------- #
library(devtools)
# install_version("ergm", version = "3.11.0", repos = "http://cran.us.r-project.org")
library(sna)
library(network)
library(igraph)
library(ergm)


setwd("C:/Users/Viv/Desktop/Lezione/AnalysisSocialNetworkLugano/Practicals/Practical3")
list.files() # Listing the files in the working directory

# -------------------------------- G(n,p) model ------------------------------ # 
set.seed(42) 
myblue <- rgb(128,177,211, max=255)

# Erdos-Renyi model
g <- sample_gnm(1000, 4950, directed = FALSE, loops = FALSE)
plot(g, vertex.color=myblue, 
     vertex.label=NA,
     vertex.size=5, 
     edge.color="black",
     layout = layout_with_mds(g))


# Gilbert model - Bernoulli random graph model
g.1 <- sample_gnp(1000, 0.05)
plot(g.1, vertex.color=myblue, 
     vertex.label=NA,
     vertex.size=5, 
     edge.color="black",
     layout = layout_with_mds(g))


# For large n the two models have the same properties
# (e.g. Poisson distribution)
par(mfrow=c(1,2))
hist(degree(g), col="lightblue",
     xlab="Degree", ylab="Frequency", main="")
hist(degree(g.1), col="lightblue",
     xlab="Degree", ylab="Frequency", main="")



# ----------------------- Preferential-attachment model ---------------------- #  
g.ba <- sample_pa(1000, directed=FALSE, out.pref=TRUE, m=1)
plot(g.ba, vertex.color=myblue, vertex.label=NA,
     vertex.size=degree(g.ba), edge.color="black")
hist(degree(g.ba), col="lightblue",
     xlab="Degree", ylab="Frequency", main="")

g.ba <- sample_pa(1000, directed=FALSE, out.pref=TRUE, m=3)
plot(g.ba, vertex.color=myblue, vertex.label=NA,
     vertex.size=degree(g.ba), edge.color="black")
hist(degree(g.ba), col="lightblue",
     xlab="Degree", ylab="Frequency", main="")


# ----------------------------- Small-world model ---------------------------- #
# Generate a regular lattice, i.e., a small-world network with rewiring prob.=0 
# The network has 50 actors, each of them with 3 neighbours 

set.seed(1908)
g.lat20 <- sample_smallworld(1, 20, 3, 0)
plot(g.lat20, layout=layout.fruchterman.reingold(g.lat20))
transitivity(g.lat20)
mean_distance(g.lat20)

g.lat20 <- sample_smallworld(1, 20, 3, 0.3)
plot(g.lat20)
transitivity(g.lat20)
mean_distance(g.lat20)

g.lat20 <- sample_smallworld(1, 20, 3, 1)
plot(g.lat20)
transitivity(g.lat20)
mean_distance(g.lat20)

par(mfrow=c(1,1))
# ----------------------------- Data description ----------------------------- #
# The folder "Practical3" contains data from a network study of 
# corporate law partnership carried out in a Northeastern US corporate law firm 
# during the period 1988-1991 in New England. The data were collected over the 
# 71 attorneys (partners and associates) of this firm. 
# Here, we consider only the subset of the 36 partners. 
# The folders includes the following files:
# - Eladv36.dat: adjacency matrix of the advice network
# - Elfriend36.dat: adjacency matrix of the friendship network.
#
# Various partners' attributes are also part of the data set. 
# The file Elattr.dat contains information on
# - id: identifier
# - office: location of the office in which the lawyer work 
#   (1=Boston; 2=Hartford; 3=Providence)
# - seniority: the years spent working for the firm
# - school: law school in which the lawyer studied 
#   (1: Harvard, Yale; 2: ucon; 3: other) 

# The full data description is available at 
# https://www.stats.ox.ac.uk/~snijders/siena/



# ------------------------------- Loading data ------------------------------- #
# - adjacency matrices 
advice <- as.matrix(read.table("ELadv36.dat",header=FALSE))
friendship <- as.matrix(read.table("ELfriend36.dat",header=FALSE))
rownames(advice) <- 1:nrow(advice)
colnames(advice) <- 1:nrow(advice)
rownames(friendship) <- 1:nrow(advice)
colnames(friendship) <- 1:nrow(advice)

# - demographic characteristics of the lawyers
attr <- read.table("ELattr36.dat",header=TRUE)
str(attr)


# ----------------------------------- CUGs ----------------------------------- #
# We test whether ties between lawyers in the same office are more likely 
# to occur in the observed network than in an Erdos-Renyi model
set.seed(1908)
homo <- function (x,attrib){
  x <- network(x)
  x %v% "school" <- attrib
  summary(x ~ nodematch("school"))
}

simHom <- cug.test(advice, homo, mode = "digraph", cmode ="edges", diag = FALSE,
                   reps = 3000, FUN.args = list(attrib=attr$school))
simHom

# EXERCISE: Compute the observed value and the empirical p-value by hand




# Similarly for reciprocity and transitivity
# Reciprocal ties are more likely to occur in the observed network 
# than in an Erdos-Renyi model
recip <- function (x) {summary(x ~ mutual)}
simRec <- cug.test(advice,recip, mode = "digraph", cmode ="edges",
                 diag = FALSE, reps = 3000)
simRec
plot(simRec)

# Transitive ties are more likely to occur in the observed network 
# than in an Erdos-Renyi model
trans <- function (x) {summary(x ~ ttriple)}
simTrans <- cug.test(advice, trans, mode = "digraph", cmode ="edges",
                 diag = FALSE, reps = 3000)
simTrans
plot(simTrans)


# EXERCISE: consider the friendship network. 
# Repeat the analysis and test for reciprocity and transitivity. 
# What can we conclude?