# Diffusion simulation models in R
# Author: Jake Fisher

# This lab will focus on how to simulate a diffusion process in R.  We will
# cover:
# - Setting up a diffusion simulation
# - Changing different parameter settings
# - Presenting the results

# Following best practices, we begin by setting up the workspace
rm(list = ls())
setwd("C:/Users/fishe/Box Sync/Home Folder jcf26/SNH2017/Instructor dropbox")
library(statnet)
library(tidyverse)
library(magrittr)
library(ggnetwork)

# Dedicated diffusion packages
library(EpiModel)
library(netdiffuseR)

# Typically, the purpose of a diffusion simulation is to illustrate the effects
# of a new diffusion process, or a new type of network structure, on the
# ultimate results of a diffusion process.  As such diffusion simulations are 
# highly variable -- after all, the point for each one is to show something new!

# However, despite the variability, most diffusion simulations build on a
# common set of building blocks:
# - an underlying network structure,
# - a transmission process, and
# - a running timer.
# In essence, we are defining an algorithm, or a set of rules, for how people
# get "infected" with a social contagion, and looking at how those rules perform
# under differnt conditions.

# To illustrate this, I will show a few examples of diffusion models.  These are
# just an introduction to what you can do; the sky is the limit!

##### Example 1: Disease epidemics ####
# The most straightforward and common use for a diffusion simulation is to
# simulate an epidemic spreading through a population.  Some examples of this
# include:
# Moody, James. 2002. "The Importance of Relationship Timing for Diffusion." 
#      Social Forces 81(1):25-56.
# Merli, M.Giovanna, James Moody, Joshua Mendelsohn, and Robin Gauthier. 2015.
#      "Sexual Mixing in Shanghai: Are Heterosexual Contact Patterns Compatible 
#      With an HIV/AIDS Epidemic?" Demography 52(3):919-42.

# In the example here, we will simulate an epidemic spreading through a random
# network.  This roughly approximates a "random mixing" scenario, where people
# interact randomly.

# To simulate random mixing, we construct the following rules:
# 1. Set up a world with a fixed number of people and a small fraction of 
#    initial infections, arranged randomly as a network
# 2. Pick an edge at random
# 3. If the edge is discordant (a susceptible-infected connection), flip a coin
#    to determine whether the non-infected person gets infected
# 4. Repeat steps 1 and 2 until everyone is infected, or until a certain number
#    of steps has passed

# First, in our simple example, let's set our parameters ahead of time
n.people <- 100
p.infection <- 0.5
pct.starting.infected <- 0.01
max.time <- 5000
contact.rate <- 0.05 / 2  # prob. of contact between any 2 people in the network

### Step 1: Set up world ###
# Create a random graph, where edges are placed randomly.  This is called a
# Bernoulli or an Erdos-Renyi random graph
set.seed(919)
net <- rgraph(n = n.people, tprob = contact.rate) %>%
  symmetrize %>%  # Make symmetric -- doesn't matter who is connected to who
  network(matrix.type = "adjacency")  # convert to statnet

# Chose a subset of people who are infected initially
infected <- sample(
  x = c(T, F),      # people can be infected (T) or susceptible (F)
  size = n.people,  # create a vector that is n.people long
  replace = T,      # with replacement, so that >1 person can be infected
  prob = c(pct.starting.infected, 1 - pct.starting.infected)
  )

### Step 2: Choose an edge ###
# For each step, we're going to choose an edge at random, and then, if the edge
# is discordant, flip a coin to determine whether the susceptible person gets
# infected.

# First, create an edgelist...
el <- as.edgelist(net) %>%
  as.data.frame %>%
  set_names(value = c("from", "to")) %>%
  
  # ... attach the values of infected...
  mutate(from.infected = infected[from], to.infected = infected[to],
         
         # ... and create a discordant variable
         discordant = (from.infected != to.infected))

# Next, choose an edge at random
random.edge <- sample(nrow(el), size = 1)

# Check if the edge is discordant
el[random.edge, "discordant"]  # it's not, so we do nothing.

# For the example, I'm going to speed this up by choosing a discordant edge
discordant.edge <- sample(which(el$discordant), size = 1)

### Step 3: Flip a coin to see if the person gets infected ###

# Now, flip a coin to see if the uninfected person gets infected
el[discordant.edge, ]  # Person 62 is the suceptible, but we will want to be
                       # able to determine that without looking manually

# A little tricky indexing to pull out the ID of the susceptible person...
who.susceptible <- with(
  el[discordant.edge, ],
  c(from, to)[!c(from.infected, to.infected)]
  )

# Flip the coin to determine if infection spreads (it does)
(infected[who.susceptible] <- sample(c(T, F), size = 1, 
                                    prob = c(p.infection, 1 - p.infection)))

### Step 4: Repeat ###
# To repeat this process, we actually embed steps 1 and 2 in a loop.

# Set up a list with the output
infections <- vector(length = max.time, mode = "list")

# Save what we already did as the first iteration
infections[[1]] <- infected

# Quick aside -- what did we create?
head(infections)

# Drop the "from.infected", "to.infected", and "discordant" columns from el,
# because they'll actually change with every iteration
el %<>% select(-from.infected, -to.infected, -discordant)

# Next, run the loop
set.seed(27708)
for (t in 2:max.time) {
  infections[[t]] <- infections[[t - 1]]
  
  # Pick an edge at random
  random.edge <- sample(nrow(el), size = 1)
  
  # If the edge is discordant, flip a coin to decide if the infection spreads
  if (with(el[random.edge, ], 
           infections[[t]][from] != infections[[t]][to])) {
    
    who.susceptible <- with(
      el[random.edge, ],
      c(from, to)[!c(infections[[t]][from], infections[[t]][to])]
      )
    
    infections[[t]][who.susceptible] <- sample(
      c(T, F), 
      size = 1, 
      prob = c(p.infection, 1 - p.infection)
    )
  }
}


# Now we have a list of who was infected at what time point.  Let's combine
# that into a data.frame, so we can work with it more easily
(results <- infections %>%
  lapply(FUN = as.data.frame) %>%
  lapply(FUN = set_names, value = "infected") %>%
  bind_rows(.id = "t") %>%
  mutate(id = rep(1:network.size(net), times = max.time),
         t = as.integer(t)) %>%
  tbl_df)

# This dataset is the raw results of our simulation, but it's usually easier to
# look at a summary.  Let's look at the number of people infected over time
infections.by.time <- results %>%
  group_by(t) %>%
  summarize(n.infections = sum(infected)) %>%
  mutate_each(funs(as.numeric), t, n.infections)

# Aside: in R, there are many ways to solve the same problem.  We could have
# gone directly from the raw data to the summaries using the apply functions:
# infections.by.time <- data.frame(
#   t = 1:max.time,
#   n.infected = sapply(infections, sum)
# ) 

# Plotting this relationship shows the 
qplot(data = infections.by.time, x = t, y = n.infections, geom = "line")

# Or, alternatively, we could look at whether people who are more central
# are infected sooner
time.to.infection <- results %>%
  group_by(id) %>%
  summarize(time.infected = min(t[infected])) %>%
  arrange(id) %>%
  mutate(indegc = degree(net))

ggplot(time.to.infection, aes(x = indegc, y = time.infected)) + 
  geom_point() + 
  geom_smooth(method = "lm")

# For fun, let's plot a few slices of the network in time
set.seed(330)
net.layout <- ggnetwork(net) %>% 
  rename(id = vertex.names)

net.layout.by.time <- split(results, f = results$t) %>%
  lapply(FUN = right_join, y = net.layout, by = "id") %>% 
  bind_rows

net.layout.by.time %>% 
  filter(t %in% c(1, 100, 250, 500, 750, 1000)) %>%
  ggplot(aes(xend = xend, yend = yend, x = x, y = y)) + 
  geom_edges(color = "lightgray") +
  geom_nodes(aes(color = infected)) + 
  facet_wrap(~ t) + 
  theme_blank()

# We could animate this with gganimate... but it doesn't work on R v. 3.4.  See
# https://raw.githubusercontent.com/jalapic/rmeetup_examples/master/ggnetwork_gganimate.R
# for example code.

# The statnet team developed a new package, called EpiModel, to simplify running
# epidemic diffusion models.  Tutorials are online here:
# http://www.epimodel.org/, and we'll walk through one to illustrate.

# Set the parameters ahead of time
param <- param.icm(
  inf.prob = 0.2,  # probability of contracting an infection in a step
  act.rate = 0.25  # speed with which people interact
  )
init <- init.icm(
  s.num = 500,  # number of people who are infected initially
  i.num = 1     # number of people who are susceptible initially
  )
control <- control.icm(
  type = "SI",  # type of model, vs., e.g., SIR or SIS
  nsims = 10,   # number of simulations to run
  nsteps = 300  # number of steps per simulation
  )

# Run the model.  ICM means individual contact model (vs. a deterministic model,
# where the results are determined by a set of differential equations)
set.seed(919)
mod <- icm(param, init, control)

# Now we can look at some summaries of the model
# Look at the state of the model at a particular time slice
summary(mod, at = 125)  

# Look at a dataset with the state of the model at a given time
mod %>%
  as.data.frame(mod = "mean") %>%
  tbl_df

# Get a plot of the number of infections by time
plot(mod)

# We can customize the plot, showing the individual simulations, for example
plot(mod, y = "i.num", sim.lines = TRUE, mean.smooth = FALSE, 
     qnts.smooth = FALSE)

# The package is quite flexible.  It can incorporate demographic changes,
# different network structures, and even custom functions for each of the 
# components of the model (like different mortality schedules).  Check the 
# website for more examples of how to do this.

##### Example 2: Diffusion of innovations #####
# The second example is the diffusion of innovations.  The process by which
# innovations spread is less well-understood.  We usually simulate it with a 
# process that's very similar to diseases diffusing.

# In this case, we will focus on a "complex contagion", meaning something you
# have to hear about from multiple people before you will adopt it.

# Setup
# 1. Construct a small world network
# 2. Choose one person at random, and set him/her and his/her neighbors as 
#    initial adopters
# 3. Infect all the people who have more than tau neighbors infected
# 4. Set the people who were infected at the previous round as "recovered"
# 5. Repeat steps 

### Set parameters in advance ###
n <- 50
tau <- 2 / 6
max.time <- 50

### Step 1: Construct network ###
# Construct a small world network
set.seed(919)
sw.net <- igraph::sample_smallworld(1, n, 2, p = .2) %>%
  intergraph::asNetwork(.) 

### Step 2: Choose a person & neighbors at random as an initial adopter ###
# Set up vector to indicate adoption
adopters <- rep(F, n)

# Choose a person at random
initial.adopter <- sample(seq_len(n), size = 1)

# Get the list of people they're attached to
initial.neighbors <- get.neighborhood(sw.net, initial.adopter)

# Set them all as "adopters"
adopters[c(initial.adopter, initial.neighbors)] <- T

### Step 3: infect all people who have more than tau neighbors infected ###

# Let's look at one person, person 18, to see how this works
ego.extract(sw.net, ego = 20, neighborhood = "out")

# Person 18 hasn't adopted
adopters[20]

# About half of their neighbors have adopted
adopters[c(18, 19, 22, 25)]
mean(adopters[c(18, 19, 22, 25)])

# To decide whether the person adopts, we test whether the fraction of
# adopters is greater than tau
mean(adopters[c(18, 19, 22, 25)]) >= tau  # adoption!

# We can update everyone simultaneously using matrix multiplication
adj.mat <- sw.net[, ]
diag(adj.mat) <- 0  # set the diagonal to 0, b/c people don't weight themselves

# We could take the sum...
adj.mat %*% adopters

#... or the percentage
(adj.mat.rn <- adj.mat / rowSums(adj.mat))
adj.mat.rn %*% adopters

# And then we calculate the people who are above our threshold
(adj.mat.rn %*% adopters) >= tau

# Note that this actually allows people to abandon the innovation if enough of
# of their neighbors are not adopters.  For now we don't want that to happen, 
# so we'll only test people who are not yet adopters
ifelse(adopters, TRUE, ((adj.mat.rn %*% adopters) >= tau))

### Step 4: Repeat ###
# Again, we can take care of this by wrapping it in a loop
adopt <- vector(mode = "list", length = max.time)
adopt[[1]] <- adopters

for (t in 2:max.time) {
  adopt[[t]] <- ifelse(adopters, TRUE, ((adj.mat.rn %*% adopt[[t - 1]]) >= tau))
}

# Note that again we get the characteristic S-shaped curve:
data.frame(
  t = 1:max.time,
  n.adopt = sapply(adopt, sum)
) %>%
  ggplot(aes(x = t, y = n.adopt)) + 
  geom_line()

# Let's plot a few frames

set.seed(330)
sw.net.layout <- ggnetwork(sw.net) %>% 
  rename(id = vertex.names)

sw.net.layout.by.time <- adopt %>%
  lapply(FUN = as.data.frame) %>% 
  lapply(FUN = set_names, value = "adopter") %>% 
  lapply(FUN = mutate, id = 1:n) %>%
  lapply(FUN = right_join, y = sw.net.layout, by = "id") %>% 
  bind_rows(.id = "t") %>% 
  mutate(t = as.integer(t))

sw.net.layout.by.time %>% 
  filter(t < 10) %>% 
  ggplot(aes(xend = xend, yend = yend, x = x, y = y)) + 
  geom_edges(color = "lightgray") +
  geom_nodes(aes(color = adopter)) + 
  facet_wrap(~ t) + 
  theme_blank()

# There's a new package for analyzing the diffusion of innovations.  It focuses
# mainly on creating survival models for diffusion data (and has a different
# data structure) so we won't cover it here.  But it's well worth checking out
# the vignette:
vignette("analyzing-medical-innovation-data")


##### Example 3: Averaging ideas #####
# A third common approach use for diffusion simulation is to suggest how
# people's beliefs change through interaction with others.  A longstanding
# assumption is that people take the weighted average of their friends.

# Set the parameters in advance
self.weight <- 0.2
max.time <- 5

# Load the data and save a few values for convenience
load("faux_add_health_with_faux_attitude.Rdata")

n.people <- network.size(add.health)  
attitude <- (add.health %v% "faux.attitude")
adj.mat <- add.health[, ]  # save the adjacency matrix to its own object

# Plotting the network shows that attitude values are randomly distributed in
# this network
set.seed(27708)
net.layout <- ggnetwork(add.health)

p <- ggplot(net.layout, aes(x = x, y = y, xend = xend, yend = yend)) + 
  geom_edges(color = "lightgray", 
             arrow = arrow(length = unit(6, "pt"), type = "closed")) + 
  theme_blank() + 
  scale_color_distiller(palette = "Spectral")

p + geom_nodes(aes(color = faux.attitude), size = 4)

# In general, we assume that people update their attitude using the average of
# their neighbors' opinions, and their own opinion.  For example, person 3 names
# 4 friends: people 6, 29, 34, and 40
ego.extract(add.health, ego = 3, neighborhood = "out")

# Those people have 4, 2, 3, and 2 as their attitudes
attitude[c(6, 29, 34, 40)]

# Their mean value is 2.75
mean(attitude[c(6, 29, 34, 40)])

# Then we assume that people take the weighted average of their own attitudes
# and their friends' attitudes
self.weight * attitude[3] + (1 - self.weight) * mean(attitude[c(6, 29, 34, 40)])

# We could repeat this for everyone individually, but this can actually be done
# by matrix multiplication, which makes it much easier to program.

# First, we have to row-normalize (set rows so that they sum to 1 by dividing
# each row by its sum) the adjacency matrix
diag(adj.mat) <- 0  # set the diagonal to 0

(adj.mat.rn <- adj.mat / rowSums(adj.mat))  # note the missing values!
                                            # Isolates => dividing by 0

# NaN does not play nicely with other values, so swap it for 0 in the row-
# normalized matrix.
adj.mat.rn[which(is.nan(adj.mat.rn), arr.ind = T)] <- 0

# In general, we assume that isolates do not change under this model.
# Mechanically, we can either remove them from the analysis or we can set their
# self-weights to 1.

# Identify outdegree isolates
isos <- (degree(add.health, cmode == "outdegree") == 0)

# Set self-weights to 1:
self.weight.vec <- rep(self.weight, times = n.people)
self.weight.vec[isos] <- 1

# Remove isolates
self.weight.no.iso <- self.weight.vec[!isos]
adj.mat.no.iso <- adj.mat.rn[!isos, !isos]
attitude.no.iso <- attitude[!isos]

# Now run the simulation...

# ... with self-weights of isolates = 1
I <- diag(n.people)  # identity matrix
SW <- diag(self.weight.vec)  # diagonal matrix with self-weights on diagonal

(adj.mat.rn %*% (I - SW) + SW) %*% attitude

# ... removing the isolates
I.no.iso <- diag(sum(!isos))
SW.no.iso <- diag(self.weight.vec[!isos])

(adj.mat.no.iso %*% (I.no.iso - SW.no.iso) + SW.no.iso) %*% attitude.no.iso

# Similarly to the other simulations, we can repeatedly perform this procedure
# using a for loop:
sim.attitude <- vector(mode = "list", length = max.time)
sim.attitude[[1]] <- attitude

for (t in 2:max.time) {
  sim.attitude[[t]] <- (adj.mat.rn %*% (I - SW) + SW) %*% sim.attitude[[t - 1]]
}

# Now let's plot the values by time, to show how the diffusion process unfolds

# Convert each list element to a data frame that ggplot can use, & combine them
layout.with.attitudes <- sim.attitude %>%
  lapply(FUN = as.data.frame) %>%
  lapply(FUN = set_names, value = "attitude.value") %>%
  lapply(FUN = mutate, vertex.names = seq_len(n.people)) %>%
  lapply(FUN = right_join, y = net.layout, by = "vertex.names") %>%
  bind_rows(.id = "iteration")

# Plot it, with facets to show how this changes over each iteration
p %+% 
  layout.with.attitudes +
  geom_nodes(aes(color = attitude.value), size = 4) +
  facet_wrap(~ iteration)
  
# A package that's under development, latentnetDiffusion, takes care of the
# matrix multiplication and reshaping of the data for you
devtools::install_github("jcfisher/latentnetDiffusion")
library(latentnetDiffusion)
library(latentnet)

# For example, to run the same simulation as above, use:
(sim.attitudes.2 <- degroot(W = add.health, Y = as.matrix(attitude), 
                            all.iter = T, self.weight = 0.8) %>%
  tidyDegroot(id = seq_len(n.people)) )

# Or, we can fit a latent space model, and simulate the diffusion process over
# the latent space:
# fit <- latentnet::ergmm(add.health ~ euclidean(d = 2))  # takes a minute!
# save(fit, file = "ergmm_2d_fit.Rdata")
load("ergmm_2d_fit.Rdata")

(ls.attitudes <- ergmmDegroot(fit, Y = as.matrix(attitude), simulate = F, 
                             all.iter = T) %>%
  tidyDegrootList(id = seq_len(n.people)))


# Plot it, to show the difference
layout.with.ls.attitudes <- ls.attitudes %>%
  rename(vertex.names = id) %>% 
  {split(., f = .$iteration)} %>% 
  lapply(FUN = right_join, y = net.layout, by = "vertex.names") %>%
  bind_rows(.id = "iteration")

p %+% 
  layout.with.ls.attitudes +
  geom_nodes(aes(color = pred.value), size = 4) +
  facet_wrap(~ iteration)
  
# The same model, simulated:
set.seed(919)
(ls.attitudes.sim <- ergmmDegroot(fit, Y = as.matrix(attitude), simulate = T, 
                                  all.iter = T) %>%
    tidyDegrootList(id = seq_len(n.people)))
  
# Plotted again
layout.with.ls.attitudes.sim <- ls.attitudes.sim %>%
  rename(vertex.names = id) %>% 
  {split(., f = .$iteration)} %>% 
  lapply(FUN = right_join, y = net.layout, by = "vertex.names") %>%
  bind_rows(.id = "iteration")

p %+% 
  layout.with.ls.attitudes.sim +
  geom_nodes(aes(color = pred.value), size = 4) +
  facet_wrap(~ iteration)