# Randomization and generating random graphs
# Author: Jake Fisher
# Set up workspace
rm(list = ls())
setwd("C:/Users/fishe/Dropbox/SN&H Files/r_scripts")
library(statnet)
library(dplyr)
# Often in statistics we are interested in a reference distribution that our
# sample, or case, is drawn from. For networks, those reference distributions
# are usually given by a distribution of random graphs. We will explore some
# simple random graph distributions today, and tomorrow Jeff will discuss more
# complicated random graph distributions with ERGMs.
# This lab will cover:
# 1. Generating random graphs
# 2. QAP tests
# 3. Running latent space models with latentnet
##### Generating random graphs #####
### Erdos-Renyi random graphs ###
# The simplest random graph distribution is the Erdos-Renyi random graph. We
# construct that by flipping a coin for every edge. statnet has a command for
# generating random graphs:
(er1 <- rgraph(
n = 50, # 50 people in the graph
m = 1, # generate 1 graph -- if greater than 1, it returns a list of matrices
tprob = .05 # probability of a tie (density of the graph)
))
# This is not magic -- we could have done the same thing by generating a matrix
# of random 1's and 0's with a given probability:
rbinom(
size = 1, # 1 binomial trial -- a Bernoulli trial
n = 50 ^ 2, # number of 0's or 1's to generate
p = .05
) %>%
matrix(., nrow = 50, ncol = 50) # structure 0's and 1's as a matrix
# Often, we'll use different random graphs to recreate a t-statistic of sorts --
# we want to see how unlikely the statistic that we observed is under a given
# null model. Let's try that for the reciprocity of the Add Health network.
# Load the Add Health data:
load("add_health_cmty1.Rdata") # saved R file from yesterday
# Store the density and network size, which we will use for the probability of
# a tie and the size of the in the random graph, respectively.
(prob <- gden(add.health))
(n <- network.size(add.health))
# Store the reciprocity of the network, which we will use for comparison later
(ah.recip <- grecip(add.health, measure = "dyadic.nonnull"))
# Generate 100 random graphs, convert them to network objects, and calculate
# the reciprocity for each
er.draws <- rgraph(n = 72, m = 100, tprob = prob) %>%
apply(., 1, network, matrix.type = "adjacency") %>%
sapply(., grecip, measure = "dyadic.nonnull")
# (As an aside, note apply loops over an index of an array. rgraph creates a
# 100 x 72 x 72 matrix, so we loop over the first index.)
# Now plot a histogram of the reciprocities of the random graph, and add a line
# showing how the observed value compares:
plot(hist(er.draws), xlim = c(0, .5))
abline(v = ah.recip, col = "red") # *way* higher than the random graphs!
# Aside: ggplot version:
# library(ggplot2)
# data.frame(recip = er.draws) %>%
# qplot(data = ., x = recip, geom = "histogram") +
# geom_vline(xintercept = ah.recip, color = "red")
#
# note the + in ggplot is very similar to %>% in dplyr. In an interview, Hadley
# said that if he had known about the %>% when he was developing ggplot2, he
# wouldn't have used the +
### Random graphs from a mixing matrix ###
# We rarely believe that a completely random graph, with no structure at all, is
# the right comparison. (If we did, studying networks wouldn't tell us much!)
# So, instead, we usually use a more sophisticated null model.
# As a next step, we consider models where people mix randomly in proportion to
# how much they mix in the observed data. The rates at which different groups
# mix are given by the mixing matrix.
(mix.mat <- mixingmatrix(add.health, "grade"))
# We can use this to get the probability that a given tie will go from one
# person to a person of another grade -- it's just the row percentages
mix.mat$matrix / rowSums(mix.mat$matrix)
(tie.probs <- prop.table(mix.mat$matrix, margin = 1)) # equivalent
# And now, we have to draw ties randomly with that probability. We'll make a
# matrix of probabilities, and then will flip a coin (draw a Bernoulli trial)
# for each possible tie.
# First, construct an empty matrix:
prob.matrix <- matrix(NA, ncol = network.size(add.health),
nrow = network.size(add.health))
# Next, fill in the values based on the mixing matrix:
grades <- add.health %v% "grade"
grade.values <- sort(unique(grades))
# The general idea is to do this...
prob.matrix[which(grades == 7), which(grades == 7)] <- tie.probs[1, 1]
# ...but that's going to require 36 lines of code. Let's write a loop instead.
for (i in 1:nrow(tie.probs)) {
for (j in 1:ncol(tie.probs)) {
prob.matrix[which(grades == grade.values[i]),
which(grades == grade.values[j])] <- tie.probs[i, j]
}
}
# And now we can flip a coin for each element of that matrix:
rg <- rgraph(n = network.size(add.health), tprob = prob.matrix)
# Ordinarily, we would do this procedure repeatedly, calculate the network
# statistic, and compare the randomly generated network statistics with the
# one from the raw data. Let's try that again with reciprocity.
mix.recip <- rgraph(n = network.size(add.health), m = 100,
tprob = prob.matrix) %>%
# Convert them all to a list of network objects
apply(., 1, network, matrix.type = "adjacency") %>%
# calculate the reciprocity for each one
sapply(., grecip, measure = "dyadic.nonnull")
# And now plot it against the true value:
plot(hist(mix.recip), xlim = c(0, .5))
abline(v = ah.recip, col = "red")
### Random graphs from a fixed degree distribution ###
# As a next step, we consider null models where the degree distribution matches
# the degree distribution in the observed data. Fortunately, this is
# implemented in igraph:
# Calculate the degree distribution from the original network:
(indegree <- degree(add.health, cmode = "indegree"))
# We actually want a table indicating the count of how many times each degree is
# observed in the real data. We use table() to do that.
(indegree <- table(indegree)) # 4 isolates, 4 people with one nomination, etc.
# Same thing for outdegree:
outdegree <- table(degree(add.health, cmode = "outdegree"))
# To generate random directed graphs with a given indegree and outdegree
# sequence, there have to be the same number of indegree values as outdegree
# values, and they must sum to the same thing. The latter happens necessarily,
# but the former does not. So we pad the smaller of the two with 0 values.
length(indegree) # 13
length(outdegree) # 11
outdegree[12:13] <- 0 # pad the smaller list with 0's
# Now we pass both of those to the function in igraph
igraph::sample_degseq(out.deg = outdegree, in.deg = indegree, method = "simple")
# Aside: converting between igraph and statnet
# igraph and statnet do very similar things, but occasionally there is something
# that you can do in one package and not in the other (like generating random
# graphs with a given degree sequence). In those cases, it's helpful to be able
# to convert from one to the other. The package to do this is called
# intergraph.
# intergraph::asNetwork() # converts to statnet from igraph
# intergraph::asIgraph() # converts to igraph from statnet
# Let's use this approach again to compare the reciprocity of the observed
# network versus the reciprocity of some simulated networks matching the same
# degree distribution:
degseqRecip <- function(indeg, outdeg) {
# Generates a random network with a given degree sequence, and returns the
# reciprocity of that network
#
# Args:
# indeg: a numeric vector of counts of people with degree 0, 1, 2...
# outdeg: a numeric vector of counts of people with degree 0, 1, 2...
# must be same length as indeg, and sum to the same value
#
# Returns:
# a numeric value, the percentage of reciprocated ties
igraph::sample_degseq(out.deg = outdegree, in.deg = indegree,
method = "simple") %>%
igraph::reciprocity(.) %>%
return(.)
}
# replicate allows you to repeat the same command repeatedly. It's useful for
# generating random numbers
more.reciprocities <- replicate(n = 100, degseqRecip(indegree, outdegree))
plot(hist(more.reciprocities), xlim = c(0, .7))
abline(v = ah.recip, col = "red")
##### QAP tests #####
# Another common null model is the permutation test. If we randomly shuffle
# the rows and columns of the adjacency matrix, how often would a correlation
# as large as the one we observed occur?
# For a simple example, let's look at whether there is a significant correlation
# between smoking behavior in the faux Add Health network and friendship.
# To begin, we construct a square matrix where the (i, j) cell is a dummy
# variable indicating whether i and j have the same smoking behavior (i.e.,
# they both smoke, or the they both don't smoke.) To construct that, we use the
# outer command (like an outer product)
smoking <- outer(
add.health %v% "PSMOKES",
add.health %v% "PSMOKES",
FUN = "=="
)
# We compare that with the matrix of friendships:
friends <- as.matrix(add.health)
# The command for QAP tests is netlm or netlogit (depending on whether you use
# a linear or logistic regression). Here we use logistic:
(qap.out <- netlogit(y = smoking, x = friends, nullhyp = "qap", reps = 100))
# You can add additional covariates by entering them as an array. Let's
# add same-grade:
grade <- outer(
add.health %v% "grade",
add.health %v% "grade",
FUN = "=="
)
# Make a 3 dimensional array -- the abind command from the abind package also
# does this.
X <- array(, dim = c(2, nrow(grade), ncol(grade)))
X[1, , ] <- friends
X[2, , ] <- grade
(qap2 <- netlogit(y = smoking, x = X, nullhyp = "qap", reps = 100))
##### Optional: latent space models #####
# Finally, another way to control for the underlying network structure is the
# latent space model. Simply put, the idea behind a latent space model is that
# we place pepole in an underlying latent "space", and then treat the ties
# between them as independent, given their positions in the space. (You can
# think of the space as the positions in a force-directed layout in a network
# picture.) This allows us to test for homophily effects, for example, or to
# generate random networks that match an observed network with very little
# complexity in model fitting.
# The most user-friendly latent space model software is in the latentnet package
# (more recent models are provided by the amen package).
library(latentnet)
latent.fit <- ergmm(add.health ~ euclidean(d = 2))
# Can add additional dimensions...
latent.fit.3d <- ergmm(add.health ~ euclidean(d = 3))
# ... latent groups ...
latent.fit.3d <- ergmm(add.health ~ euclidean(d = 2, G = 5))
# ... or homophily effects
latent.fit.3d <- ergmm(add.health ~ nodematch("grade") + euclidean(d = 2))
