Introduction to Software: R

• Flexibility
  
• Data
• Convert formats
• Customized cleaning
• New methods
• Slow GUI development
• Open source community
• Can write your own functions
• Replication
• Personal record keeping
• Collaboration
• Transparency
• Awareness

• Open source
• Community
  
• Large and growing
  
• Made by statisticians
• Wealth of packages (i.e., libraries)
• Features
• Documentation
• Moderately easy language
• Graphics

• Shared platform
• Networks
• Statistics
• Data management
• Graphics
• Active network software development
• Network metrics
  
• igraph
  
• sna
• tnet
  
• Network modeling
• statnet
• RSiena

• Command line-only interface
• Reliable, Replicable, Records
• GUI options do exist, but...
• Outside of scope for our purposes
• Assignment. Functionally equivalent methods

• x <- 1 # Preference: primary method 

• x = 1 # Preference: within functions    

• 1 -> x # Preference: never    

• Once assigned, x becomes an "object"

• ls()    

Cheburashka <- "a cartoon character"Cheburashka
cheburashka

• R runs within a folder on your hard drive
  
• Location where R reads & writes
  
• Which folder are you in now?

orig.folder <- getwd()orig.folder
dir.create("new.folder")
list.files()
setwd("new.folder")list.files()
setwd(orig.folder)
file.remove("new.folder")list.files()

• Preloaded (e.g., base, stats)
• Dependencies
• User-contributed
• Social network analysis
  
• statnet
• igraph
• tnet
• RSiena
• ...
• Graphics
• Advanced regression models
  
• Web scraping
• Spatial
• ...

 install.packages("igraph")

 library(igraph)

 detach("package:igraph", unload=TRUE)

x <- NA

NA # Non-applicable, missing
NULL # Undefined or empty value
NaN # Not a Number
Inf # Infinite 

is.na(x)
is.null(x)
is.nan(x)
is.infinite(x)

• Elements
  
• Numeric
• Character
• Use double quotes
• Logical
• TRUE, FALSE , T , F , 1 ,  0
• Objects
  
• Vectors
• Matrices
• Arrays
• Data frames
• Lists
• Closures
  

is.numeric(42)
is.numeric("42")
is.numeric(3.12)
is.logical(T)
is.logical(FALSE)
is.character(FALSE)
is.character("42")
is.character("HSE Network Workshop") 

as.numeric("42") + 1
as.character(as.numeric("0") + 1)
as.character(!as.logical("1"))
as.character(as.logical(as.numeric("0") + 1))
as.character(is.logical(as.logical(as.numeric("0")) + 1))
as.character(!as.logical(as.numeric("1")))  

RU.cit <- c("Mos", "StP", "Nov", "Yek") # Create vector
RU.cit <- c(RU.cit, "Niz") # Add element
RU.cit[3] # Retrieve element 3
RU.pops <- c(10.4, 4.6, 1.4, 1.3) # Create vector
names(RU.pops) <- RU.cit[c(1, 2, 3, 4)] # Name vector
RU.pops[c(4,2,3,1)] # Permute the order and retrieve
length(RU.pops) # Retrieve the length
RU.vis <- c(TRUE, TRUE, FALSE, FALSE) #Create vector
RU.vis[-c(4, 2)] # Remove elements 4 and 2

all(1:5 == c(1, 2, 3, 4, 5))
all(1:5 == seq(from = 1, to = 5, by = 1)) 

rep(1, times = 5)
rep(NA, times = 4)
rep("Beetlejuice", times = 3)

• Two dimensions
  
• nrow(); ncol()  
• Assignment
  
• x[1, 5] <- 1; x[2, ] <- x[1, 5] > 0  
• More columns, rows
  
• cbind(); rbind()
• Names:  rownames();  colnames()

• Three dimensions
• dim(); length(x[1, 1, ])  
• Assignment
  
• x[3, 2, 6] <- 1; x[4, , 2] <- NA  
• Names:  dimnames()

RU.mat <- cbind(pop = RU.pops, vis = RU.vis)
rownames(RU.mat) <- RU.cit[c(1:4)]
colnames(RU.mat)

RU.mat <- matrix( , nrow = length(RU.pops), ncol = 2)
colnames(RU.mat) <- c("pop", "vis")
rownames(RU.mat) <- names(RU.pops)
RU.mat[ , "pop"] <- RU.pops
RU.mat[ , "vis"] <- RU.vis 

RU.mat[1, ]RU.mat[ , "pop"]
RU.mat[1, 2]RU.mat[1, , ] # What does the error mean?

library(compiler)InvHypSineTrans.cmp <- cmpfun(InvHypSineTrans)system.time(InvHypSineTrans(0:10000000))system.time(InvHypSineTrans.cmp(0:10000000))

1. Find the randomized normal distribution function
   
2. Use it to create two vectors of length 100
   
1. One with mean of zero & standard deviation of one
1. Assign it x
2. A second with mean of 50 & standard deviation of 25
1. Assign it y
3. Find the linear regression model command
   
4. Regress y on x

test.val <- "It works!"
test.val2 <- paste(test.val, "Really")

source("test.val.source.txt")
ls()
print(test.val)
print(test.val2)

• History
  
• RHistory files retain a log of your session
• May be opened in a text editor
• Should end in ".RHistory" extention
• Usage

  savehistory("Filename.RHistory")
  loadhistory("Filename.RHistory") 

• Images
  
• Saves every object in the environment (i.e., ls())
• Usage

• save.image("Filename.RData")
  load("Filename.RData")
  

• dput(), dget()
• Saves and retrieves individual objects

• Network science tradition
• General and complex network metrics
• Data storage
• Speed and efficiency

• General social network metrics
• Data storage
• Modeling and statistical inference

• Many identical function names
• Different data formats

• asIgraph()
• asNetwork()

• Weighted
• Two-mode
• Longitudinal

• Longitudinal network modeling
• Longitudinal behavioral network modeling
• Cross-sectional network modeling
• p* / ERGM

• Matrix
• Edge list
• Pajek

• Actor dataframe
• Edge dataframe

library(igraph)

wikidat <- read.table(file = "http://pastebin.com/raw.php?i=UVmTznBj", sep = "\t", header = TRUE, skip = 6)
# 'sep' is the separator
# 'header' indicates column names
# 'skip' are the number of initial rows to ignore

wikidat <- graph.data.frame(d = wikidat, directed = TRUE)
summary(wikidat) # Nodes: 7115; Edges: 103689

library(igraph)

marvnet <- read.graph("http://bioinfo.uib.es/~joemiro/marvel/porgat.txt", format = "pajek")
summary(marvnet)

• Case
• Early study on electronic information exchange
• Precursor to email, late-1970s
• Subjects
• 32 researchers
• Cross-discipline social network analysts
• Relationships, expressed in matrices
• Heard of - Met - Friend - Close friend at two time points
• Message frequency between time points
• Attributes
• Vertex: discipline
• Vertex: citation count
• Edges: affinity
• Edges: message frequency
  

firstdir <- getwd()
newdir <- paste(firstdir, "EIES_data", sep = "/")
download.file("https://www.stats.ox.ac.uk/~snijders/siena/EIES_data.zip", "EIES_data.zip", method = "wget")
unzip("EIES_data.zip", exdir = newdir)
setwd(newdir)

eiesatts <- read.table("eies.dat", as.is = TRUE, stringsAsFactors = FALSE, strip.white = TRUE, skip = 3, quote = "", nrows = 32)
colnames(eiesatts) <- c("Researcher", "Original_Code", "Citations", "Discipline")
eiesatts$Discipline_Name <- c("Sociology", "Anthropology", "Statistics", "Psychology")[eiesatts$Discipline]
eiesatts$Discipline_Name[3] <- "Communication"
eiesatts$Discipline_Name[27] <- "Mathematics"

friend1 <- read.table("EIES.1", sep = " ", as.is = TRUE, strip.white = TRUE)
friend2 <- read.table("EIES.2", sep = " ", as.is = TRUE, strip.white = TRUE)
eiesmessages <- read.table("messages.asc", nrows = 32, strip.white = TRUE)
eies <- list(t1 = friend1, t2 = friend2, messages = eiesmessages)
rm(friend1, friend2, eiesmessages)
eies <- lapply(eies, as.matrix)

eiesvatts <- function(g, atts = eiesatts){
    V(g)$name <- as.character(1:vcount(g))
    for (i in 1:ncol(atts))
        g <- set.vertex.attribute(g, name = colnames(atts)[i], value = atts[, i])
    return(g)
    }
eies <- lapply(eies, eiesvatts)
rm(eiesvatts, eiesatts)

sapply(list.files(), file.remove)
setwd(firstdir)
file.remove(newdir)
file.remove("EIES_data.zip")
rm(firstdir, newdir)

?read.graph
?graoh.data.frame
?graph.adjacency
?graph.adjlist
?graph.constructors
# Query database for popular data sets
?nexus.list
?nexus.search
?nexus.get
?nexus.info

is.simple(wikidat) # Verify it's simple
is.directed(wikidat) # Verify it's directed
vcount(wikidat) # Number of vertices
ecount(wikidat) # Number of edges
graph.density(wikidat) # Density

V(wikidat)$indegree <- degree(wikidat, mode = "in")
V(wikidat)$outdegree <- degree(wikidat, mode = "out")
V(wikidat)$totaldegree <- degree(wikidat, mode = "total")

all.vatts <- list.vertex.attributes(wikidat)[-1]
all.vatts <- sapply(all.vatts, get.vertex.attribute, graph = wikidat)
summary(all.vatts)
hist(log(V(wikidat)$outdegree), main = "Wiki Admin Votes", xlab = "Outdegree (log)")

eies.closefriends.t1 <- subgraph.edges(eies$t1, eids = which(E(eies$t1)$etype == "Close Friend"), delete.vertices = FALSE)
eies.closefriends.t1 <- as.undirected(eies.closefriends.t1)
eies.closefriends.t2 <- subgraph.edges(eies$t2, eids = which(E(eies$t2)$etype == "Close Friend"), delete.vertices = FALSE)
eies.closefriends.t2 <- as.undirected(eies.closefriends.t2)
eies.closefriends <- list(t1 = eies.closefriends.t1, t2 = eies.closefriends.t2)
rm(eies.closefriends.t1, eies.closefriends.t2)

remove.isolates <- function(g){
    nonisolates <- which(degree(g) > 0)
    return(induced.subgraph(g, nonisolates))
    }
eies.closefriends.nonisolates <- lapply(eies.closefriends, remove.isolates)
rm(remove.isolates)

unlist(dyad.census(wikidat))
reciprocity(wikidat)

• Edgewise?
• Dyadic?
• Dyadic, non-null ("ratio")?

transitivity(eies.closefriends$t1) # Unidrected network
V(eies.closefriends$t1)$loc.trans <- transitivity(eies.closefriends$t1, "local")
triad.census(wikidat) # Directed network
source("http://pastebin.com/raw.php?i=J4fPjFuN")
transitivity_directed(wikidat) # Not in igraph

V(eies.closefriends$t1)$name[which.max(V(eies.closefriends$t1)$loc.trans)]

average.path.length(wikidat)
diameter(wikidat)
V(wikidat)$name[(get.diameter(wikidat))]
hist(shortest.paths(wikidat), main = "Histogram of Shortest Path Lengths", xlab = "Geodesic Lengths")

V(wikidat)$betw <- betweenness(wikidat)E(wikidat)$eb <- edge.betweenness(wikidat)hist(log(E(wikidat)$eb), main = "Histogram of Edge Betweenness", sub = "Wiki Admin Voting", xlab = "Edge Betweenness (log)")

#How many weak and strong components do we have? 
sapply(c("weak", "strong"), no.clusters, graph = wikidat)

#Notice the distributions
clusters(wikidat, mode = "weak")$csize
clusters(wikidat, mode = "strong")$csize

#Assign membership as attributesV(wikidat)$comp.w <- clusters(wikidat, mode = "weak")$membership
V(wikidat)$comp.s <- clusters(wikidat, mode = "strong")$membership# Are the vertexes part of the giant component?
V(wikidat)$giantcomp.w <- V(wikidat)$comp.w == which.max(clusters(wikidat, mode = "weak")$csize)
V(wikidat)$giantcomp.s <- V(wikidat)$comp.s == which.max(clusters(wikidat, mode = "strong")$csize)

V(wikidat)$closeness <- closeness(wikidat)
V(wikidat)$evcent <- evcent(wikidat)$vector

V(wikidat)$kc.undir <- graph.coreness(as.undirected(wikidat, mode = "collapse"))
#How are undirected k-cores related to centrality?
kc.cent.corr.fun<-function(x, y = V(wikidat)$kc.undir){
  a <- get.vertex.attribute(wikidat, name = x)
  b <- y
  cor.est <- cor.test(a, b, method = "spearman", exact = FALSE)$estimate
  return(cor.est)
  }
kc.cent.corr.fun <- cmpfun(kc.cent.corr.fun)

wikidat.comms <- multilevel.community(as.undirected(wikidat, mode = "collapse"))
wikidat.comms$modularity
V(wikidat)$comms <- wikidat.comms$membership
xtabs(~V(wikidat)$comms) # 3, 4, 5, 9, and 11 are the biggest

wikidat.comms.w <- multilevel.community(as.undirected(wikidat, mode = "collapse"), weights = max(E(wikidat)$eb) - E(wikidat)$eb)
wikidat.comms.w$modularity # Slightly higher
V(wikidat)$comms.w <- wikidat.comms.w$membership
xtabs(~V(wikidat)$comms.w) # 4, 5, 14, and 15 are the biggest

save.image("NetSchool.RData")
savehistory("NetSchool.RHistory")

plot(eies.closefriends$t1)

library(RColorBrewer)library(compiler)
my_ecdf_fun <- function(x.val, all.vals){
    ret.val <- mean(x.val >= all.vals)
    return(ret.val)
    }
my_ecdf_fun <- cmpfun(my_ecdf_fun)
vsize <- sapply(V(eies.closefriends$t1)$Citations, my_ecdf_fun, all.vals = V(eies.closefriends$t1)$Citations)
vsize <- 7.5 + 15 * vsize

each.discipline <- unique(V(eies.closefriends$t1)$Discipline_Name)
n.disciplines <- length(each.discipline)
vcol <- brewer.pal(n.disciplines, name = "Accent")
vcol <- vcol[match(V(eies.closefriends$t1)$Discipline_Name, each.discipline)]

plot(eies.closefriends$t1, vertex.color = vcol, vertex.size = vsize, vertex.label.color = "black", vertex.label.family = "sans")

vcount(wikidat) # 7115 nodes
ecount(wikidat) # 103689 edges
plot(wikidat)

sampsize <- round(vcount(wikidat) * 0.05) # 5% sample
subsamp <- sample.int(vcount(wikidat), sampsize)
wikidat.sg <- induced.subgraph(wikidat, subsamp)
plot(wikidat.sg)

sampsize <- round(vcount(wikidat) * 0.05) # 5% sample
subsamp <- sample.int(vcount(wikidat), sampsize)
wikidat.sg <- induced.subgraph(wikidat, subsamp)
plot(wikidat.sg)

• Nodes too large
• Labels don't provide information

plot(wikidat.sg, vertex.size = 3, vertex.label = "")

subsamp <- which(degree(wikidat.sg, mode = "all") > 0)
wikidat.sg <- induced.subgraph(wikidat.sg, subsamp)
plot(wikidat.sg, vertex.size = 3, vertex.label = "", edge.arrow.size = 0.5)

• Grossly underrepresents degree
• Majority of alters absent
• Underestimates connectivity
• Conceptually naïve
• Sample could reflect structure
• Replication limitations

# Sample based on betweenness and large communities
cond1 <- V(wikidat)$betw > 1
cond2 <- V(wikidat)$comms %in% c(3, 4, 5, 9, 11)
subsamp <- which(cond1 & cond2)
wikidat.sg <- induced.subgraph(wikidat, subsamp)
vcount(wikidat.sg) # From 7,115 to 1,372 nodes
ecount(wikidat.sg) # From 103,689 to 40,587 edges
# Still too big to see anything.
# Systematic 25% sample on eigenvector centrality
evrank <- rank(V(wikidat.sg)$evcent, "random")
nsyssamp <- seq(from = 1, to = max(evrank), by = 4)
nsyssamp <- which(evrank%in% nsyssamp)
wikidat.sg <- induced.subgraph(wikidat.sg, nsyssamp)
# Systematic 25% sample on edge betweenness
ebrank <- rank(E(wikidat.sg)$eb, "random")
esyssamp <- seq(from = 1, to = max(ebrank), by = 4)
esyssamp <- which(ebrank %in% esyssamp)
wikidat.sg <- subgraph.edges(wikidat.sg, esyssamp)
plot(wikidat.sg, vertex.size = 3, vertex.label = "", edge.arrow.size = 0.25)

?layout.auto
wikidat.sg.lo <- layout.fruchterman.reingold(wikidat.sg, params = list(niter = 2500, area = vcount(wikidat.sg)^3))
V(wikidat.sg)$x <- wikidat.sg.lo[,1]
V(wikidat.sg)$y <- wikidat.sg.lo[,2]
plot(wikidat.sg, vertex.size = 3, vertex.label = "", edge.arrow.size = 0.25)

• Sizes
  
• Colors and shading
• Categorical variables
• Interval and ordinal variables
• Transparency

• Nodes
• Non-continuous: Shape, node color, border color
• Continuous: Node size, border width, color
• Edges
• Non-continuous: Color, line type, arrowhead type
• Continuous: Width, color
• Labels: Size, color, visibility

ewidth <- sapply(E(wikidat.sg)$eb, my_ecdf_fun, E(wikidat.sg)$eb)
E(wikidat.sg)$width <- (1.25 * ewidth)^2
plot(wikidat.sg, vertex.size = 3, vertex.label = "", edge.arrow.size = 0.25, edge.curved = TRUE)

vcols <- brewer.pal(3, "Accent")
names(vcols) <- unique(V(wikidat.sg)$comms)
V(wikidat.sg)$color <- vcols[as.character(V(wikidat.sg)$comms)]
V(wikidat.sg)$size <- 6 * sapply(V(wikidat.sg)$evcent, my_ecdf_fun, V(wikidat.sg)$evcent)
plot(wikidat.sg, vertex.label = "", edge.arrow.size = 0.25, edge.curved = TRUE)

ebrank <- rank(E(wikidat.sg)$eb, "random")
esyssamp <- seq(from = 1, to = max(ebrank), by = 2)
esyssamp <- which(ebrank %in% esyssamp)
wikidat.sg <- subgraph.edges(wikidat.sg, esyssamp)
plot(wikidat.sg, vertex.label = "", edge.arrow.size = 0.25, edge.curved = .25)

par(mfrow = c(2, 2), mar = c(0.1, 0.1, 0.1, 0.1))
plot(wikidat.sg, vertex.label = "", edge.arrow.size = 0.25, edge.curved = .25, layout = layout.kamada.kawai, margin = c(0, 0, 0, 0))
plot(wikidat.sg, vertex.label = "", edge.arrow.size = 0.25, edge.curved = .25, layout = layout.circle, margin = c(0, 0, 0, 0))
plot(wikidat.sg, vertex.label = "", edge.arrow.size = 0.25, edge.curved = .25, layout = layout.svd, margin = c(0, 0, 0, 0))
plot(wikidat.sg, vertex.label = "", edge.arrow.size = 0.25, edge.curved = .25, layout = layout.mds, margin = c(0, 0, 0, 0))

wikidat.sg.elist <- get.edgelist(wikidat.sg, names = FALSE)
edgecolor <- apply(wikidat.sg.elist, 1, function(x){
    c1 <-V(wikidat.sg)$comms[x[1]]
    c2 <- V(wikidat.sg)$comms[x[2]]
    return(ifelse(c1 == c2, vcols[as.character(c1)], "darkgrey"))
    })
E(wikidat.sg)$color <- edgecolordev.off()
plot(wikidat.sg, vertex.label = "", edge.arrow.size = 0.25, edge.curved = .25, layout = layout.svd, margin = c(0, 0, 0, 0))

• Less is more
• Fewer modifications are better than more
• Smaller data is easier to represent than larger data
• Transparency when reporting visualization
• Informative and intuitive over beautiful
• Practical over gimmicky
  
• Science over art

?igraph.plotting

source("http://pastebin.com/raw.php?i=SErkr5sg")

• ei_index
  
• transitivity_directed
  
• uman_game
  
• rewire_dyads
  

graph.order <- seq(from = 10, to = 100, by = 5)
graph.trans <- sapply(graph.order, function(x) replicate(1000, transitivity(erdos.renyi.game(x, .1))))
graph.apl <- sapply(graph.order, function(x) replicate(1000, average.path.length(erdos.renyi.game(x, .1))))
par(mfrow = c(1, 2))
plot(x = graph.order, y = apply(graph.trans, 2, mean), ylim = c(0, 1), xlim = c(0, 100), type = "l", xlab = "Order", ylab = "Transitivity", main = "Tie Prob = 0.10; Simulations=1000")
lines(x = graph.order, y = apply(graph.trans, 2, quantile, probs = .025, na.rm = TRUE), col = "red", lty = 2, type = "l")
lines(x = graph.order, y = apply(graph.trans, 2, quantile, probs = .975, na.rm = TRUE), col = "red", lty = 2, type = "l")
plot(x = graph.order, y = apply(graph.apl, 2, mean), ylim = c(0, 5), xlim = c(0, 100), type = "l", xlab = "Order", ylab = "Average Path Length", main = "Tie Prob = 0.10; Simulations = 1000")
lines(x = graph.order, y = apply(graph.apl, 2, quantile, probs = .025, na.rm = TRUE), col = "red", lty = 2, type = "l")
lines(x = graph.order, y = apply(graph.apl, 2, quantile, probs = .975, na.rm = TRUE), col = "red", lty = 2, type = "l")

graph.tieprob <- seq(from = .01, to = .30, by = .01)graph.trans <- sapply(graph.tieprob, function(x) replicate(1000, transitivity(erdos.renyi.game(100, x))))graph.apl <- sapply(graph.tieprob, function(x) replicate(1000, average.path.length(erdos.renyi.game(100, x))))par(mfrow = c(1, 2))plot(x = graph.tieprob, y = apply(graph.trans, 2, mean), ylim = c(0, 1), xlim = c(0, .3), type = "l", xlab = "Tie Probability", ylab = "Transitivity", main = "Order = 100; Simulations = 1000")lines(x = graph.tieprob, y = apply(graph.trans, 2, quantile, probs = .025, na.rm = TRUE), col = "red", lty = 2, type = "l")lines(x = graph.tieprob, y = apply(graph.trans, 2, quantile, probs = .975, na.rm = TRUE), col = "red", lty = 2, type = "l")plot(x = graph.tieprob, y = apply(graph.apl, 2, mean), ylim = c(0,8), xlim = c(0,.3), type = "l", xlab = "Tie Probability", ylab = "Average Path Length", main = "Order=100; Simulations = 1000")lines(x = graph.tieprob, y = apply(graph.apl, 2, quantile, probs = .025, na.rm = TRUE), col = "red", lty = 2, type = "l")lines(x = graph.tieprob, y = apply(graph.apl, 2, quantile, probs = .975, na.rm = TRUE), col = "red", lty = 2, type = "l")

• Reassignment methods
• Graph generators

• Permutation tests
• Rewiring tests

• Permute an attribute to test its relationship to a nodal property
• Permute a graph to test its relationship to another graph

# Do conservative blogs have more links than liberal ones?
polblogs <- nexus.get(id = "polblogs") #0 = liberal, 1 = conservative

polblogs.totdeg <- degree(polblogs, mode = "all")
obs.diff <- mean(polblogs.totdeg[V(polblogs)$LeftRight == 1]) - mean(polblogs.totdeg[V(polblogs)$LeftRight == 0]) 

perm.fun <- function(){
  perm.attribs <- sample(V(polblogs)$LeftRight)
  return(mean(polblogs.totdeg[perm.attribs == 1]) - mean(polblogs.totdeg[perm.attribs == 0]))
  }; perm.fun<-cmpfun(perm.fun)
exp.diffs <- replicate(1000, perm.fun())
mean(obs.diff > exp.diffs) # ~0.86, not significantly greater than by chance assignment

obs.val <- ei_index(polblogs, V(polblogs)$LeftRight)
exp.val <- replicate(1000, ei_index(polblogs, sample(V(polblogs)$LeftRight)))
mean(obs.val < exp.val) # We expect fewer outgroup ties

# Friendship homophily by discipline among researchers
# Time 1
obs.val.t1 <- ei_index(eies.closefriends$t1, V(eies.closefriends$t1)$Discipline_Name)
exp.val.t1 <- replicate(1000, ei_index(eies.closefriends$t1, sample(V(eies.closefriends$t1)$Discipline_Name)))

# Time 2
obs.val.t2 <- ei_index(eies.closefriends$t2, V(eies.closefriends$t2)$Discipline_Name)
exp.val.t2 <- replicate(1000, ei_index(eies.closefriends$t2, sample(V(eies.closefriends$t2)$Discipline_Name)))
# Compare the homophily effect & its non-random probability
c(t1 = obs.val.t1, t1.p = mean(obs.val.t1 < exp.val.t1), t2 = obs.val.t2, mean(obs.val.t2 < exp.val.t2))

• Shuffle edges (or dyads)
• Shuffle edge weights

karate <- nexus.get("karate")
karate.lo <- layout.kamada.kawai(karate, params = list(niter = 1000))

par(mfrow = c(1, 2))
plot(karate, vertex.size = 5, vertex.label = V(karate)$name, vertex.label.family = "sans", margin = c(0, 0, 0, 0), edge.curved = .33, layout = karate.lo)

karate.rw <- rewire(karate, mode = "simple", niter = 10*ecount(karate))
plot(karate.rw, vertex.size = 5, vertex.label = V(karate)$name, vertex.label.family = "sans", margin = c(0, 0, 0, 0), edge.curved = .33, layout = karate.lo)

obs.trans <- transitivity(karate)karate.rw.fun <- function(){  krw <- rewire(karate, mode = "simple", niter = 10 * ecount(karate))  return(transitivity(krw))  }; karate.rw.fun <- cmpfun(karate.rw.fun)exp.trans <- replicate(1000, karate.rw.fun())mean(obs.trans > exp.trans) # ~0.87

n.researchers <- vcount(eies$t1)
n.rewires <- 10 * n.researchers * (n.researchers - 1)
obs.val.t1 <- transitivity_directed(eies$t1)
obs.val.t2 <- transitivity_directed(eies$t2)

# Will take some time.
# Would recommend at least 1000, not 250, replications
exp.val.t1 <- replicate(250, transitivity_directed(rewire_dyads(eies$t1, n.rewires)))
exp.val.t2 <- replicate(250, transitivity_directed(rewire_dyads(eies$t2, n.rewires)))

c(t1 = obs.val.t1, t1.p = mean(obs.val.t1 > exp.val.t1), t2 = obs.val.t2, t2.p = mean(obs.val.t2 > exp.val.t2))

1. Erdős–Rényi model
   
2. Power Law model
   
3. Watts and Strogatz model
4. Exponential random graph models
5. ...and many others

• Exact #ties or density
• Directed, undirected
• Dyad census-conditioned

wikidat.in <- degree(wikidat, mode = "in")
wikidat.out <- degree(wikidat, mode = "out")
wikidat.plfit.in <- power.law.fit(wikidat.in)
wikidat.plfit.out <- power.law.fit(wikidat.out)
print(c(wikidat.plfit.in$alpha, wikidat.plfit.out$alpha))
exp.incent.pl.fun <- function(){
  new.net <- static.power.law.game(vcount(wikidat), ecount(wikidat), exponent.out = wikidat.plfit.out$alpha, exponent.in = wikidat.plfit.in$alpha)
  return(centralization.degree(new.net, mode="in", loops = FALSE)$centralization)
  }; exp.incent.pl.fun <- cmpfun(exp.incent.pl.fun)
exp.incent.pl <- replicate(500, exp.incent.pl.fun())
mean(obs.incent > exp.incent.pl) # Obs more centralized
par(mfrow = c(1, 2))
hist(exp.incent.ec, main = "E-R Indegree Centralization")
hist(exp.incent.pl, main = "Power Law Indegree Centralization") # ~10x more centralized graphs

• Begin with a lattice
  
• Each actor initially has the same degree, neighborhood index (nei) * 2
• Rewire with a probability of p  
  

#The following two plots should be approximately equal
dev.off(); par(mfrow = c(1, 2))
sw.net <- watts.strogatz.game(dim = 1, size = 100, nei = 1, p = .1)
plot(g1.rewire, vertex.size = 5, vertex.label = "", main = "Rewired Lattice")
plot(sw.net, vertex.size = 5, vertex.label = "", layout = g1.layout, main = "Watts-Strogatz")

par(mfrow = c(1, 2))
g2.ws <- watts.strogatz.game(dim = 1, size = 100, nei = 1, p = .9)
g2.er <- erdos.renyi.game(100, (2 * 100 / 2) / (100 * 99 / 2), directed = FALSE)
plot(g2.ws, vertex.size = 5, vertex.label = "", main = "Watts-Strogatz")
plot(g2.er, vertex.size = 5, vertex.label = "", main = "Erdős–Rényi")

• Clustering follows a different mechanism
• Distance can take on a different meaning
• Depends upon conceptualization of tie strength

• Theoretical problems
• Assumes an "event" creates ties among all participating actors
• Do you always know everyone at each party you attend?
• Some actors and events are much more inclusive
• Should their relationships hold the same weight?
• Methodological problems
• Creates a complete subgraph ("clique") around each event
• Increases density

• Limited functions in igraph (R, Python, C, Ruby)
  
• More capabilities in tnet (R)
  
• NetworkX (Python)
  

library(igraph)
library(tnet)
library(compiler)
marvnet <- read.graph("http://bioinfo.uib.es/~joemiro/marvel/porgat.txt", format = "pajek")
V(marvnet)$name <- V(marvnet)$id

marvnet.tnet <- get.edgelist(marvnet, names = FALSE)
marvnet.tnet <- as.tnet(marvnet.tnet, type = "binary two-mode tnet")

• Limited functions in igraph (R, Python, C, Ruby)
  
• More capabilities in tnet (R)  and NetworkX (Python)  

download.file("https://raw.githubusercontent.com/kjhealy/revere/master/data/PaulRevereAppD.csv", "PaulRevere.csv", method = "wget")
amrev <- read.csv("PaulRevere.csv", row.names = 1)
file.remove("PaulRevere.csv")
amrev.el <- which(amrev == 1, arr.ind = TRUE)
amrev.el[, 1] <- rownames(amrev)[as.numeric(amrev.el[, 1])]
amrev.el[, 2] <- colnames(amrev)[as.numeric(amrev.el[, 2])]
amrev <- graph.data.frame(amrev.el, directed = FALSE)
V(amrev)$type <- V(amrev)$name %in% amrev.el[, 2]

amrev.tnet <- get.edgelist(amrev, names = FALSE)
amrev.tnet <- as.tnet(amrev.tnet, type = "binary two-mode tnet")

ra_clust <- function(g){
    g.motif <- graph.motifs(g, size = 4)
    threepaths <- g.motif[7]
    fourcycles <- 4 * g.motif[9]
    ret.val <- fourcycles / (fourcycles + threepaths)
    return(ret.val)
    }
ra_clust <- cmpfun(ra_clust)

amrev.raclust <- ra_clust(amrev)
amrev.raclust

# Much slower, but faster version available (email Opsahl)
amrev.oclust <- clustering_tm(amrev.tnet)
amrev.oclust

local_ra_clust <- function(g, vind){
    new.g <- graph.neighborhood(g, order = 2, nodes = vind)[[1]]
    new.g.motif <- graph.motifs(new.g, size = 4)
    threepaths <- new.g.motif[7]
    fourcycles <- 4 * new.g.motif[9]
    retval <- ifelse(fourcycles + threepaths > 0, fourcycles / (fourcycles + threepaths), NaN)     return(retval)
    }
local_ra_clust <- cmpfun(local_ra_clust)
amrev.raclust.loc <- sapply(1:vcount(amrev), local_ra_clust, g = amrev)

# Identify both the actors and those actors with valid clustering coefficients
amrev.raclust.loc.actors.ind <- which(is.nan(amrev.raclust.loc[which(V(amrev)$type == FALSE)]) == FALSE)
amrev.raclust.loc.actors <- data.frame(name = V(amrev)$name[amrev.raclust.loc.actors.ind])
amrev.raclust.loc.actors$clust <- amrev.raclust.loc[amrev.raclust.loc.actors.ind]
amrev.raclust.loc.actors$rank <- rank(amrev.raclust.loc.actors$clust)
amrev.raclust.loc.actors[order(amrev.raclust.loc.actors$clust, amrev.raclust.loc.actors$name), ]

• Random network generation
  
• Edge rewires

# Generate a random network with the same parameters as the Marvel network
n.primary <- sum(V(marvnet)$type == FALSE)
n.secondary <- vcount(marvnet) - n.primary
marvnet.rand.el <- rg_tm(ni = n.primary, np = n.secondary, ties = ecount(marvnet))
head(marvnet.rand.el) # It's an edgelist.
# Export it to igraph
summary(marvnet.rand.el)
# Because both the first and second mode begin at 1 ...
marvnet.rand.el[, "p"] <- marvnet.rand.el[, "p"] + n.primary
marvnet.rand.type <- c(rep(FALSE, times = n.primary), rep(TRUE, times = n.secondary))
marvnet.rand.atts <- data.frame(name = 1:vcount(marvnet), type = marvnet.rand.type)
marvnet.rand <- graph.data.frame(marvnet.rand.el, directed = FALSE, vertices = marvnet.rand.atts)

par(mfrow = c(2, 2))
tempdegree <- degree(marvnet, v = which(V(marvnet)$type == FALSE))
hist(log(1 + tempdegree), main = "Observed", xlab = "Issues in Which Each Character Appears (log)", ylab = "Frequency")
tempdegree <- degree(marvnet.rand, v = which(V(marvnet)$type == FALSE))
hist(log(1 + tempdegree), main = "Random", xlab = "Issues in Which Each Character Appears (log)", ylab = "Frequency")

tempdegree <- degree(marvnet, v = which(V(marvnet)$type == TRUE))
hist(log(1 + tempdegree), main = "Observed", xlab = "Characters Appearing in Each Issue (log)", ylab = "Frequency")
tempdegree <- degree(marvnet.rand, v = which(V(marvnet)$type == TRUE))
hist(log(1 + tempdegree), main = "Random", xlab = "Characters Appearing in Each Issue (log)", ylab = "Frequency")

• Good for small, dense networks
  
• Highly general
• Controls for size and density
• Does not control for degree distribution

# tnet method.
# Do not use now, as it will take A VERY LONG TIME.
# Also allows for weighted two-mode data and can shuffle weights

amrev.rand.rw <- rg_reshuffling_tm(amrev.tnet) # Reshuffles each edge

# igraph method. Inadviseable for large networks# First, recreate the bipartite network, but treat it as directed to preserves the bipartite structure
# Second, rewire the directed two-mode network
# Third, remove the direction
rewire_tm <- function(g, niter = NULL){
    g.el <- get.edgelist(g)
    new.g <- graph.data.frame(g.el, directed = TRUE)
    V(new.g)$type <- V(g)$type
    nswaps <- ifelse(is.null(niter), 10 * ecount(g), niter)
    new.g <- rewire(new.g, mode = "simple", nswaps)
    new.g <- as.undirected(new.g)
    return(new.g)
    }
rewire_tm <- cmpfun(rewire_tm)

1. Project observed two-mode data to one-mode
2. Take graph-level measurement
3. Simulate
1. Simulate two-mode network from observed two-mode networks
2. Project simulated two-mode to a one-mode network
3. Take graph-level measurement on simulated one-mode network
4. Repeat to create distribution
4. Compare observed measurement to simulated distribution

• Dichotomous (Brieger. 1974. Social Forces.)
• Sum (Brieger. 1974. Social Forces.)
• Weighted
• k threshold cut-off
• Newman. 2001. Phys Rev E.
• Collaboration weight penalized by number of collaborators
• Backbone
• Convert vertex pairs to z-scores using frequency counts
  
• Set alpha cut-off
• Retain edges beyond cut-off (+/-)
• Ahn,  Ahnert, Bagrow, & Barabási. 2011. Scientific Reports.
  
• Serrano,  Boguñá, & Vespignani. 2009 PNAS.
• Demonstration: http://pastebin.com/UNi0ryG3

# tnet
marv.dichot <- projecting_tm(marvnet.tnet, method = "binary") # Takes me ~1.5 minutes
head(marv.dichot, 15) # Performed on "primary" mode only
# igraph.
# Set 'which' to "true" for the other mode
# Set it to "both" for a list of both projections
marv.dichot <- bipartite.projection(marvnet, multiplicity = FALSE, which = "false")

• Dichotomous
• "Any commonality"
• Sum
• Weighted
• "Number of commonalities"
• k threshold cut-off
• "Any commonality beyond k"
• Newman (2001)
• "Probability of dyadic interaction"
• Backbone
• "Dyadic interaction far above (or below) 'chance'"
• "Chance" based upon the level of an individual's activity

• Dichotomous
• Simple: number of collaborations away
• Weighted (including Newman)
• Distance between adjacent vertices is 1 / edge weight
• Average path length may be less than 1!
• Normalization (Opsahl, Agneessens, and Skvoretz 2010)
• Divide the weights by the mean of all weights
• Then calculate distance

# Dichotomous
shortest.paths(marv.newman, v = 1, to = 2, weights = rep(1, ecount(marv.newman))) # 3

# Summed weights
shortest.paths(marv.newman, v = 1, to = 2, weights = E(marv.newman)$weight) # 4
shortest.paths(marv.newman, 1, 2, weight = E(marv.newman)$weight / mean(E(marv.newman)$weight)) # 1.18

# Newmanedges2keep <- which(is.infinite(E(marv.newman)$newman.w) == FALSE)
marv.newman <- subgraph.edges(marv.newman, edges2keep, delete.vertices = FALSE)
shortest.paths(marv.newman, v = 1, to = 2, weights = E(marv.newman)$newman.w) # 0.38
shortest.paths(marv.newman, v = 1, to = 2, weights = E(marv.newman)$newman.w / mean(E(marv.newman)$newman.w)) # 1.33

• Symbols
• Primary mode: circular nodes
• Secondary modes: square nodes
• Layout
• Bipartite
• Primary mode to the left
• Secondary mode to the right
• Multidimensional layouts appropriate also
• Projected as one-mode
• Any of the conventional layouts

• Character co-appearance
• Project to one-mode
• Subsample
• Popular characters (two-mode degree)
• Strong relationships (Newman's weight)
• Remove isolates
• Clean the names
• Color by cohesion
• Scale vertices by closeness
• Scale edge widths by Newman's weight

issues48 <- which((degree(marvnet) >= 48) & (V(marvnet)$type == FALSE))
issues48 <- which(V(marv.newman)$name %in% V(marvnet)$name[issues48])
marv.newman.sg <- induced.subgraph(marv.newman, issues48)
marv.newman.sg <- subgraph.edges(marv.newman.sg, which(E(marv.newman.sg)$newman.w > 4), delete.vertices = TRUE)

library(RColorBrewer)
# Community detection
V(marv.newman.sg)$comm <- multilevel.community(marv.newman.sg, weights = E(marv.newman.sg)$weight)$membership
max(V(marv.newman.sg)$comm) # 25 is too many communities to portray
# k-cores
V(marv.newman.sg)$kc <- graph.coreness(marv.newman.sg)
kcmax <- max(V(marv.newman.sg)$kc)
kcmax # 10 on a continuous scale is good.
min(V(marv.newman.sg)$kc) # Because there are no isolates, minimum is 1.
corecolors <- brewer.pal(max(V(marv.newman.sg)$kc), "Reds")
corecolors <- c(rep("white", times = kcmax - 9), corecolors)
V(marv.newman.sg)$color <- corecolors[V(marv.newman.sg)$kc]

V(marv.newman.sg)$closeness <- closeness(marv.newman.sg, weights = E(marv.newman.sg)$weight / mean(E(marv.newman.sg)$weight), normalize = FALSE)
V(marv.newman.sg)$size <- ecdf(V(marv.newman.sg)$closeness)(V(marv.newman.sg)$closeness)
V(marv.newman.sg)$size <- V(marv.newman.sg)$size * 6
V(marv.newman.sg)$label.cex <- V(marv.newman.sg)$size / (4.5 * 2)