SAOM_ABM

Aims

We are going to use RSiena to simulate a stochastic actor-oriented model as an agent-based model.

A more extensive script for simulating with RSiena is given in NetworkSimulation.R. That script also illustrates how to simulate multiple waves, somethign that might be useful if you want to track the trajectory over time.

Model assumptions

In a time-window \(t \in [0,1)\), \(n\) actors are assumed to make changes to their behaviour and/or network. Once a decision to make a change, or to consider the opportunity to change, is made by an actor, the actor makes a decision based on the current state of the system, picking the course of action that maximises a stochastic utility function.

Rates

For each behaviour \(r \in \mathcal{R}=\{1,\ldots, R\}\), each actor has a rate of change \(\lambda_{i,r}(X(t))\), where \(X(t)\) is the current state of the system at time \(t\). The time until actor \(i \in V\) gets an opportuntiy to make a change to behavior \(r\) is exponentially distributed with rate \(\lambda_{i,r}(X(t))\). This means that we can pick the actor and behaviour that might cahnge in a number of ways. One way is to draw \(n \times R\) variables \(\langle S_{i,r} \rangle_{i \in V,r \in \mathcal{R}}\) from exponential distributions with rates \(\langle \lambda_{i,r}(X(t)) \rangle_{i \in V,r \in \mathcal{R}}\), and then pick the actor \(i\) and behaviour \(r\) that has the smallest value \(S_{i,r}\), i.e. the shortest wating time. Another way of picking actor \(i\) and behaviour \(r\), is by using the property of the exponential distribution that the probability that \(S_{i,r}\) is minimum, is \[\frac{\lambda_{i,r}}{\sum_i \sum_r \lambda_{i,r}} .\] Once the actor \(i\) and behaviour \(r\) have been chosen, and a potential change has happened, time is incremented \(t := t+S_{i,r}\).

Decision process

Given that actor \(i\) is given the opportunity to make a change to the behaviour \(r\), the actor chooses the make the change that maximises \[ U_i(r,\Delta_{i,r,k}X(t))=f_i(r,\Delta_{i,r,k}X(t))+\epsilon_{i,r,k,t} \] where

  • \(\Delta_{i,r,k}X(t)\) is the state that would result from \(i\) making change \(k\) to behaviour \(r\),
  • \(U_i(r,X)\) is the utility of state \(X\) to actor \(i\) with respect to behaviour \(r\)
  • \(f_i(r,X)\) is the objective function of the utility of state \(X\) to actor \(i\) with respect to behaviour \(r\)
  • \(\epsilon_{i,r,k,t}\) is an extreme value type one variable

Typically, for a network change, \(\Delta_{i,r,k}X(t)\) is the network that would result from \(i\) making a change to their tie to \(k\) by toggling the current state. This means that \(\Delta_{i,r,k}X(t)\) would be an adjacency matrix that is identidal to that of \(X(t)\), with the exception that you set the cell \((i,j)\) to \(1-X(t)_{i,j}\). The universe of changes to the network avaible to \(i\) is \((\Delta_{i,k}: k=1,\ldots,i-1,i+1,\ldots,n)\) as well as \(\Delta_{i,i}\) which is defined as a non-change.

A simple example of a diffusion and network model

First load the RSiena, RSiena, andRSiena` packages:

library(sna)
library(network)
library(RSiena)
library(ggplot2)

We will use a network of \(n=32\) actors as a starting network (but you could use any network you like)

tmp4[is.na(tmp4)] <- 0 # remove missing
mynet1 <- sienaDependent(array(c(tmp4, tmp4), dim=c(32, 32,2)),allowOnly = FALSE)# note that we want to simulate both creation and deletion of ties

Create a binary starting behaviour

mybeh.1 <- matrix(runif(2*32)<.2,32,2)+0 # remove missing
mybeh <- sienaDependent( mybeh.1, type = "behavior" ,allowOnly = TRUE)# we only want to model an absorbing state

Create the data object

mydata <- sienaDataCreate(mynet1,mybeh)

Simulate model

Decide what effects to include and set parameter values

myeffects <- getEffects( mydata )
myeffects <- setEffect( myeffects, name = "mynet1",density,initialValue = -1.4163)# set the density paramter
##   effectName          include fix   test  initialValue parm
## 1 outdegree (density) TRUE    FALSE FALSE    -1.4163   0
myeffects <- setEffect( myeffects,recip,initialValue = 1.1383 )# set the reciprocity paramter
##   effectName  include fix   test  initialValue parm
## 1 reciprocity TRUE    FALSE FALSE     1.1383   0
myeffects <- includeEffects( myeffects, transTrip, transRecTrip )# adding two transitivity effects
##   effectName                  include fix   test  initialValue parm
## 1 transitive triplets         TRUE    FALSE FALSE          0   0   
## 2 transitive recipr. triplets TRUE    FALSE FALSE          0   0
myeffects <- setEffect( myeffects, egoX, 
                         interaction1 = "mybeh" ,initialValue = -.5 )# social selection effects
##   effectName include fix   test  initialValue parm
## 1 mybeh ego  TRUE    FALSE FALSE       -0.5   0
myeffects <- setEffect( myeffects,  altX, 
                         interaction1 = "mybeh" ,initialValue = .25 )# social selection effects
##   effectName  include fix   test  initialValue parm
## 1 mybeh alter TRUE    FALSE FALSE       0.25   0
myeffects <- setEffect( myeffects, egoXaltX,
                         interaction1 = "mybeh",initialValue = .5 )# social selection homophilye
##   effectName              include fix   test  initialValue parm
## 1 mybeh ego x mybeh alter TRUE    FALSE FALSE        0.5   0
myeffects <- setEffect( myeffects, avAlt, name="mybeh",
                         interaction1 = "mynet1",   initialValue = 1.2 ) # a positive influence effect
##   effectName          include fix   test  initialValue parm
## 1 mybeh average alter TRUE    FALSE FALSE        1.2   0
myeffects <- setEffect(myeffects, name = "mynet1", Rate, type="rate",
                            initialValue =  4.2525)
##   effectName                  include fix   test  initialValue parm
## 1 basic rate parameter mynet1 TRUE    FALSE FALSE     4.2525   0
myeffects <- setEffect(myeffects, name = "mybeh", Rate, type="rate",
                            initialValue = 3)
##   effectName          include fix   test  initialValue parm
## 1 rate mybeh period 1 TRUE    FALSE FALSE          3   0

Set up the simulation settings:

sim_model  <-  sienaAlgorithmCreate( 
                          projname = 'sim_model',
                          cond = FALSE, 
                          useStdInits = FALSE, nsub = 0 ,
                          n3 = 100,# run 100 simulations
                          simOnly = TRUE)   # by setting simOnly to TRUE RSiena won't estimate

Run simulations

sim_ans <- siena07( sim_model,#simulation settings
                    data = mydata,# our starting data
                    effects = myeffects,# the effects and paramter values we have set for model
                    returnDeps = TRUE,# this is to actually return the simulated networks and behaviours
                    batch=TRUE )    
## 
## Start phase 0 
## theta:  4.25 -1.42  1.14  0.00  0.00  0.25 -0.50  0.50  3.00  0.00  1.20 
## 
## Start phase 3

Extract simulated data

Selection and influence tables

Selection and influence tables plot the part of the utility function that corresponds to the selection and influence effects. Note that these tables are the theoretical effects and not based on what is actually simulated

Read in the script from a website

source('http://www.stats.ox.ac.uk/~snijders/siena/SelectionTables.r')

Calculate the social selection table:

vname <- "mybeh"
name <- "mynet1"
levls <- 0:1
vselect <- selectionTable(sim_ans, mydata, name, vname, levls)
## Dependent network mynet1 ; actor variable mybeh .
## Parameters found are
##     altX     egoX egoXaltX 
##     0.25    -0.50     0.50 
## Check that these parameter values are correct!
## The subtracted mean value of mybeh is 0.15625 .
vselect
##   ego vego valter      select
## 1   0    0      0  0.05126953
## 2   0    0      1  0.22314453
## 3   1    1      0 -0.52685547
## 4   1    1      1  0.14501953
sp <- ggplot(vselect, aes(valter, select, group=ego, colour=ego))
sp +geom_point() + geom_smooth(size=1.2, span=3) +
scale_colour_hue() +
scale_x_continuous(breaks=levls) +
theme(legend.key=element_blank())+
labs(x=paste(vname),
y=paste("selection function"),
title=paste("Effect",vname,"on",name),
colour=paste(vname)) +
theme_grey(base_size=26, base_family="") +
theme(legend.key.width=unit(1, "cm")) +
theme(plot.title=element_text(hjust=0.5))
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

Read in the Influence script from a website

source('http://www.stats.ox.ac.uk/~snijders/siena/InfluenceTables.r')

Calculate the social influence table:

name <- "mybeh"
zname <- "mynet1"
levls <- 0:1
zselect <- influenceMatrix(sim_ans, mydata,zname,  name, levls)
## Network mynet1 ; dependent behavior mybeh .
## Parameters found are
## avAlt 
##   1.2 
## 
## Levels of alter refer to average values of alter's behavior.
zselect
##             0          1
## 0  0.02929688 -0.1582031
## 1 -0.15820312  0.8542969
zselect <- influenceTable(sim_ans, mydata,zname,  name, levls)
## Network mynet1 ; dependent behavior mybeh .
## Parameters found are
## avAlt 
##   1.2 
## 
## Levels of alter refer to average values of alter's behavior.
sp <- ggplot(zselect, aes(zego, select, group=alter, colour=alter))
sp + geom_point() + geom_smooth(size=1.2, span=3) +
scale_colour_hue() +scale_x_continuous(breaks=levls) +
theme(legend.key=element_blank())+labs(x=paste(name),
                                       y=paste('evaluation function'),
                                       title=paste('Influence effect',zname,'on',name),
                                       colour=paste(name,'\nalter')) +
theme_grey(base_size=26, base_family="") +
theme(legend.key.width=unit(1, "cm")) +
theme(plot.title=element_text(hjust=0.5))
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

Extracting networks

The function networkExtraction returns the network as an edge list of class “network” according to the network package (used for package sna). Unfortunately there is a bug that I can’t get around so here we do it the hard way

#attr(sim_ans$f[['Data1']]$depvars[["mynet1"]],  "type")
simnets <- vector('list',100)
for (i in c(1:100)){

emptyNetwork <- network::network.initialize(32, 
            bipartite = NULL)
matrixNetwork <- sparseMatrixExtraction(i, sim_ans$f, sim_ans$sims, 
        1, 'Data1', "mynet1")
    sparseMatrixNetwork <- as(matrixNetwork, "dgTMatrix")
simnets[[i]] <- network::network.edgelist(cbind(sparseMatrixNetwork@i + 
            1, sparseMatrixNetwork@j +1, 1), emptyNetwork)
}
hist(gden(simnets))# chart density of simulated networks

For simulated behaviour, extracting the simulated values is straightforward (once you know or guess the structure of the ‘sims’ object; see ‘9.1 Accessing the generated networks’ in the RSiena Manual Manual)

#attr(sim_ans$f[['Data1']]$depvars[["mybeh"]],  "type")

simbeh <- matrix(NA,32,100)
for (i in c(1:100)){
simbeh[,i] <- sim_ans$sims[[i]][[1]][[2]][[1]]
# [[iteration]][[period]][[dependent variable]][[]]
}
 hist(colMeans(simbeh))# proportion of adopters

Let’s look at the number of homophilous and heterophilour ties

mat <- as.matrix.network(simnets[[1]])# convert edgelist to matrix
sum( mat[simbeh[,1]==1, simbeh[,1]==1] )# number 1-> 1
## [1] 60
sum( mat[simbeh[,1]==1, simbeh[,1]==0] )# number 1-> 0
## [1] 35
sum( mat[simbeh[,1]==0, simbeh[,1]==1] )# number 0-> 1
## [1] 44
sum( mat[simbeh[,1]==0, simbeh[,1]==0] )# number 0-> 0
## [1] 20

Calculate a function that calculates homophilous ties across simulations or just loop

hom.tab <- matrix(NA,100,4)
colnames(hom.tab) <- c('1=>1','1=>0','0=>1','0=>0')
for (i in c(1:100))
{
mat <- as.matrix.network(simnets[[i]])# convert edgelist to matrix
hom.tab[i,1] <-  sum( mat[simbeh[,i]==1, simbeh[,i]==1] )# number 1-> 1
hom.tab[i,2] <- sum( mat[simbeh[,i]==1, simbeh[,i]==0] )# number 1-> 0
hom.tab[i,3] <- sum( mat[simbeh[,i]==0, simbeh[,i]==1] )# number 0-> 1
hom.tab[i,4] <- sum( mat[simbeh[,i]==0, simbeh[,i]==0] )# number 0-> 0
}
plot(ts(hom.tab))