# setting working directory

# insert path to your working directory
# don't forget change \ to /
setwd("insert a path to your working directory here") 

install.packages("network") # installing library
library(network) # loading library
install.packages("sna")
library(sna)
install.packages("bipartite") # installing package for bipartite networks analysis
install.packages("statnet")
install.packages("PCIT") 

data() # available datasets

data(flo) # loading data

?flo # information about dataset

flo # matrix of Florentine family wedding ties 

gplot(flo, usearrows=F)


###----------matrix operations----------


# creates matrix manually

mat <- matrix(c(1,2,3,4,5,
                6,7,8,9,10,
                11,12,13,14,15,
                16,17,18,19,20,
                21,22,23,24,25), nrow=5, ncol=5) # first we put in data, afterwards we define # of rows and columns
mat

m <- matrix(1:25, ncol=5, nrow=5) # elements are sorted column-wise (default)
m

n <- matrix(1:25, ncol=5, nrow=5, byrow=T) # elements are sorted row-wise (byrow=T)
n

# matrix operations

t(n) # transposition
diag(n) # matrix diagonal

m1 <- matrix(rep(c(1,2),2), nrow=2, ncol=2)
m1
m2 <- matrix(rep(2,4), nrow=2, ncol=2)
m2

m1 + m2 # adding up matrices

m1 * m2 # multiplying matrices' elements

2*m1 # multiplying matrix by number

m1 %*% m2 # multiplying matrix by matrix


###----------bipartite network example----------


library(bipartite) # loading package for bipartite networks analysis

acm <- as.data.frame(read.csv("n_24_twomode_matrix.csv")) # reading in data
rownames(acm) <- acm$X
acm$X <- NULL # deleting default rownames
acm <- as.matrix(acm) # conversion to incidence (actors by concepts) matrix

head(acm) # first six rows of data set
View(acm) # displays data

### one-mode projections (concepts by concepts and actors by actors matrices)

concept.adj <- as.one.mode(acm, project="higher", weighted=T) # concepts cooccurence (adjacency) matrix 

actor.adj <- as.one.mode(acm, project="lower", weighted=TRUE) # actors overlap (adjacency) matrix

plot.sociomatrix(concept.adj) # shows as a sociomatrix

plot.sociomatrix(actor.adj) # shows as a sociomatrix

concept.net <- network(concept.adj, valued = T, bipartite = F, directed = F) # converts to one-mode network

actor.net <- network(actor.adj, valued = T, bipartite = F, directed = F) # converts to one-mode network

### bipartite network

bp.net <- network(acm, valued = T, bipartite = T, directed = F) # converts to bipartite network

gplot(bp.net, # data (network object)
      gmode="twomode", # two-mode network
      mode="fruchtermanreingold", # lay out
      interactive=F,  # non-interactive lay out
      label=network.vertex.names(bp.net), # node labels
      boxed.labels=F, label.pos=5, label.cex=0.75, # label demarcation, position and font size
      vertex.col=c(rep("lightgreen", dim(acm)[1]), rep("lightblue", dim(acm)[2])), # node colors for actors and concepts
      vertex.sides=c(rep(5, dim(acm)[1]), rep(20, dim(acm)[2])), # node shapes for actors and concepts
      usearrows=F, # without arrows (undirected ties)
      edge.col="darkgrey" # ties' color
      )

plotweb(acm) # plotting two-mode network with valued nodes and edges

visweb(as.matrix(acm)) # plotting two-mode sociomatrix

# plottting one-mode projections (cooccurence and overlap matrices)

plot.sociomatrix(concept.adj) # plotting concept by concept sociomatrix

gplot(concept.net, gmode="graph", mode="kamadakawai", usearrows=F, label=network.vertex.names(concept.net), 
      label.pos=5, vertex.col="lightblue", edge.col="darkgrey"
      #, edge.lwd=concept.adj # adds weighted ties (very messy)
      )

plot.sociomatrix(actor.adj) # plotting actor by actor sociomatrix

gplot(actor.net, mode="kamadakawai", label=network.vertex.names(actor.net), usearrows=F, 
      label.pos=5, edge.col="darkgrey", vertex.col="lightgreen", vertex.sides=5 
      #, edge.lwd=actor.adj # adds weighted ties (very messy)
      )


###----------creating a network----------


g <- network.initialize(5) # creates empty network with 5 nodes
plot(g) # plotting the network 

g[1,2] <- 1 # adds tie between nodes 1 and 2
plot(g)

g[,3] <- 1 # adds ties between node 3 and all remaning nodes
plot(g)

g[1,3] # shows value of a tie between nodes 1 and 3

d <- g[-1,-1] # creates a subgraph that excludes node 1
plot(network(d)) # plotting the subgraph

g[,] <- 0 # removes all ties = empty network
plot(g)


###----------data import and processing----------


library(statnet) # loading statnet package

m <- read.csv("matrix_2001.csv", sep=";", row.names=1) # loads relational data (organized as adjacency matrix), values seperated by ;, first column = row name
head(m) # first six rows
View(m) # displays data
class(m)

m <- as.matrix(m) # conversion into matrix
m

as.sociomatrix(m) # sociomatrix

# conversion into network / network object

net <- network(m,                    # data
               directed=FALSE,       # undirected ties
               ignore.eval=FALSE,    # saves values of ties
               names.eval="value")   # names attribute as "value"

class(net)
net # basic info

plot(net)


###----------descriptive statistics----------


summary(net)

## local indicators: centrality

degree(net, gmode="graph") # degree centrality
degree(net, gmode="graph", rescale=T)
deg <- degree(net, gmode="graph")/37 # standardization: /(n-1)
round(deg, 2)

closeness(net, gmode="graph") # result is already standardized: *(n-1)
clo <- closeness(net[c(-28,-30),c(-28,-30)], gmode="graph") # closeness centrality
clo

betweenness(net, gmode="graph") # betweenness centrality
betweenness(net, gmode="graph", rescale=T)
bet <- betweenness(net, gmode="graph")/((37*36)/2) # standardization: /((n-1)*(n-2)/2)
round(bet, 2)

library(PCIT) # includes function for clustering coefficient

localClusteringCoefficient(m) # local clustering coefficient

cc <- localClusteringCoefficient(m) 
cc[!is.finite(cc)] <- 0 # conversion of NaN values to 0
cc

## global indicators

gden(net, mode="graph") # density

connectedness(net) # connectivity

gcc <- sum(cc)/length(cc) # average clustering coefficient
gcc

dyad.census(net) # dyads
grecip(net) # reciprocity

triad.census(net, mode="graph") # triads 
gtrans(net,  mode="graph") # transitivity

clique.census(net, mode="graph") # cliques

ccen <- clique.census(net, mode="graph")
ccen$clique.count # number of n-cliques (2,3,4 cliques in this case)

centralization(net, degree, mode="graph", normalize=FALSE) # degree centralization
centralization(net, degree, mode="graph", normalize=TRUE) # degree centralization standardized

centralization(net, closeness,  mode="graph", normalize=TRUE) # closeness centralization

centralization(net, betweenness,  mode="graph", normalize=TRUE) # betweenness centralization

kcores(net, mode="graph", cmode="freeman", ignore.eval=T) # identification of k-cores


##----------vizualization----------


## weighted ties

val_net <- as.sociomatrix(net, "value") # save tie values under attribute called "value"
plot.sociomatrix(val_net) # weighted sociomatrix

## attributional data

node <- read.csv("nodes_2001.csv", header=T, sep=";", row.names=1, stringsAsFactors=F) # loading attributional data
View(node)
node <- node[-6,][-17,][-19,][-21,][-22,][-27,][-27,][-31,][-35,][-39,] # removing cases with missing data

cinc <- node[,3] # Composite Index of National Capability 
cinc

# plotting the network including the weighted ties and attributes

gplot(net                               # data (network object)
      ,gmode="graph"                    # undirected graph
      #, mode="fruchtermanreingold"
      #, mode="circle" 
      #, mode="spring"
      ,mode="kamadakawai"               # lay-out of the network
      ,label=network.vertex.names(net)  # nodes names
      #,vertex.cex=(cinc)*100           # node size is based on attribute value (CINC)
      #,vertex.cex=deg*2                 # node size is based on attribute value (degree centrality)
      ,vertex.col="lightblue"           # nodes' color (can be also based on some attribute value)
      ,label.pos=5                      # position of the label
      ,label.cex=0.75                   # font size
      ,boxed.labels=F                   # label without any demarcation
      ,vertex.sides=20                  # node shape (based on the number of sides - 20 approximates round shape)
      ,edge.col="darkgrey"              # ties' color
      ,edge.lwd=net%e%('value')*10)     # projects saved values of weighted ties onto unweighted network called "net"

# titles, legends etc. 

title(main="Interdependence network (European natural gas market (2001))", # main title
      sub="Data (COW NMC_v4.0, Eurostat, World Bank)",                     # subtitle
      cex.main=0.75,                                                       # font size of main title  
      font.main=2,                                                         # bold (2)
      cex.sub=0.75,                                                        # font size of subtitle 
      font.sub=3)                                                          # italics (3)
legend('bottomright',                                                      # legend location
       c("edge width: interdependence",                                    # text of the legend
         "node size: national capabilities"),                              
       cex=0.75,                                                           # font size of legend
       bty="n",                                                            # legend demarcation (n = none)
       title="",                                                           # legend title (none)
       text.font=1)                                                        # ordinary font (1)


##----------exploratory analysis----------


# matrix of Euclidean distances

dist <- sedist(net, method="euclidean", mode="graph") # calculation of Euclidean distances
colnames(dist) <- colnames(m) # copies column names from adjacency matrix to Euclidean distances matrix
rownames(dist) <- colnames(dist) # copies column names to row names
dist <- round(dist, 2) # round to 2 decimal places
dist
range(dist)

# cluster analysis

c_dist <- equiv.clust(dist, method="euclidean", plabels=colnames(dist)) # hierarchical clustering
plot(c_dist, main="Euclidean clustering", labels=colnames(dist)) # dendrogram vizualization
rect.hclust(c_dist$cluster, h=5) # cluster demarcation (h specifies value for demarcation)

## egocentric network

ego_net <- ego.extract(net, neighborhood="combined") # extraction of egocentric network
ego_net
gplot(ego_net$CZR, boxed.labels=F, label=colnames(ego_net$CZR), usearrows=F) # CZ
gplot(ego_net$NTH, boxed.labels=F, label=colnames(ego_net$NTH), usearrows=F) # NED
gplot(ego_net$RUS, boxed.labels=F, label=colnames(ego_net$RUS), usearrows=F) # RUS
gplot(ego_net$SLO, boxed.labels=F, label=colnames(ego_net$SLO), usearrows=F) # SLO