############################################################################### #### ShowBehavioralChains.R #### ############################################################################### # This is a script demonstrating the operation of the siena07 # # with the parameter returnChains=TRUE # # to get an output of all the simulated ministeps, # # and use this for behavioral dependent variables. # # Written by Tom Snijders, using the script WorkOnChains.r as a basis. # # version 21-03-2022 # ############################################################################### #### Set up a data set for simulation #### ############################################################################### # Take a basis for having a model with behavior from Rscript04SienaBehaviour.R. library(RSiena) friend.data.w1 <- s501 friend.data.w2 <- s502 friend.data.w3 <- s503 drink <- s50a smoke <- s50s # create dependent variables: friendship <- sienaDependent(array(c(friend.data.w1, friend.data.w2, friend.data.w3 ), dim = c(50, 50, 3 ) ) ) drinkingbeh <- sienaDependent(drink, type = "behavior" ) # construct data set: mydata <- sienaDataCreate(friendship, drinkingbeh ) # specify model: myeff <- getEffects(mydata ) myeff <- includeEffects(myeff, transTrip, transRecTrip) myeff <- includeEffects(myeff, egoX, altX, simX, interaction1 = "drinkingbeh" ) myeff <- includeEffects(myeff, name = "drinkingbeh", avAlt, interaction1 = "friendship" ) myeff algo <- sienaAlgorithmCreate(projname = NULL, seed=642) (ans <- siena07(algo, data = mydata, effects = myeff )) # Now we have constructed a sienaFit object to work with. ############################################################################### #### Get the simulated chains #### ############################################################################### # Let us simulate and return the chains: for n3=10 times: sim.algorithm <- sienaAlgorithmCreate(projname = NULL, cond=FALSE, nsub=0, n3=10, simOnly=TRUE, seed=123) sim.ans <- siena07(sim.algorithm, data=mydata, effects=myeff, returnChains=TRUE, returnDataFrame=TRUE) # Note that returnChains=TRUE is not supported with nbrNodes > 1. ############################################################################### #### Inspect what was created #### ############################################################################### theChain <- sim.ans$chain str(theChain, 1) # ten runs str(theChain[[1]], 1) # two periods str(theChain[[1]][[1]], 1) # many ministeps # The chain has the structure # chain[[run]][[group]][[period]][[ministep]] # For a single group, the [[group]] level is omitted. # Length of the chain for period 1 of run 1: length(theChain[[1]][[1]][[1]]) # 283 # all runs for period 1: sapply(theChain, function(x)length(x[[1]][[1]])) # all runs for period 2: sapply(theChain, function(x)length(x[[1]][[2]])) # the first ministep: theChain[[1]][[1]][[1]][[1]] # a nicer view of this: unlist(theChain[[1]][[1]][[1]][[1]]) # Now we inspect the chain for period 1 of run 1: sapply(theChain[[1]][[1]][[1]], length) # But there are some components that are always NULL sapply(theChain[[1]][[1]][[1]], function(x){length(unlist(x))}) # Note that in C++, indices for numbering start with 0. # For the meaning of the first 10 fields of a ministep: # From siena07utilities.cpp: # SET_STRING_ELT(colnames, 0, mkChar("Aspect")); network or behaviour; # SET_STRING_ELT(colnames, 1, mkChar("Var")); # SET_STRING_ELT(colnames, 2, mkChar("VarName")); # SET_STRING_ELT(colnames, 3, mkChar("Ego")); # SET_STRING_ELT(colnames, 4, mkChar("Alter")); # SET_STRING_ELT(colnames, 5, mkChar("Diff")); # difference in dependent behavior variable: 0 if nothing changes # SET_STRING_ELT(colnames, 6, mkChar("ReciRate")); # SET_STRING_ELT(colnames, 7, mkChar("LogOptionSetProb")); # SET_STRING_ELT(colnames, 8, mkChar("LogChoiceProb")); # SET_STRING_ELT(colnames, 9, mkChar("Diagonal")); # This is from C++ code; note that in C++, numbering starts at 0. # Therefore, for a one-mode network the values of Ego and Alter run # from 0 to n-1, where n is the number of actors. # If the mini-step is a network mini-step for a one-mode network, # the change made can be inferred by comparing ego and alter. # If ego and alter are the same node, then no change has occurred; # if ego and alter are different nodes, then the value for the tie # from ego to alter is toggled (1 -> 0; 0 -> 1). # For a two-mode network the values of Alter run from 0 to m, # where m is the number of nodes in the second mode; # here the value Alter=m means that no change has occurred. # Inspect: table(sapply(theChain[[1]][[1]][[1]], function(x){x[[1]]})) # Aspect table(sapply(theChain[[1]][[1]][[1]], function(x){x[[2]]})) # Dependent variable table(sapply(theChain[[1]][[1]][[1]], function(x){x[[3]]})) # Name of dependent variable table(sapply(theChain[[1]][[1]][[1]], function(x){x[[4]]})) # Ego table(sapply(theChain[[1]][[1]][[1]], function(x){x[[5]]})) # Alter table(sapply(theChain[[1]][[1]][[1]], function(x){x[[6]]})) # Difference if behavior ############################################################################### #### Transform for seeing the evolution of the behaviour #### ############################################################################### # Look at only the first run: runNumber <- 1 (chainLength <- length(theChain[[runNumber]][[1]][[1]])) # The second variable here is the behavioral variable: depVariable <- 2 # let us consider only period 1: period <- 1 # What is the variable we consider: mydata$depvars[[depVariable]] table(mydata$depvars[[depVariable]][,,period]) # To get the results of the ministeps, # we create functions that mimic how the dependent variables # are transformed by a ministep. # For networks, define the toggle function: toggle <- function(mat,i,j) { if (i != j) {mat[i,j] <- 1-mat[i,j]} mat } # For behavior, define the behChange function: behChange <- function(z, i, d) { z[i] <- z[i] + d z } # In script WorkOnChains.r, the networks are studied, # here only the behaviour. # The behaviour starts at (startBeh <- mydata$depvars[[2]][,,period]) ## [data for wave 1] # To calculate the intermediate behaviors # and the values of the desired statistics: # Transform the sequence of ministeps to a matrix: theSteps <- matrix(unlist(theChain[[runNumber]][[1]][[1]]), 11, chainLength) # This is only for the run with number runNumber, # although 10 runs were generated. # Inspect the first three steps: theSteps[,1:3] # For the meaning of the 11 rows, see above. # Extract the relevant rows: vars <- as.numeric(theSteps[2,])+1 # number of dependent variable egos <- as.numeric(theSteps[4,])+1 alters <- as.numeric(theSteps[5,])+1 difs <- as.numeric(theSteps[6,]) # Start with a dummy non-step egos <- c(0, egos) vars <- c(0, vars) difs <- c(0, difs) # Prepare for calculations of intermediate states: avz <- rep(NA, chainLength+1) zlist <- list() # Start calculations: # Here we start with the initial value. theBeh <- startBeh avz[1] <- mean(startBeh) zlist <- c(zlist, list(startBeh)) # Pursue calculations: for (ministep in 1:chainLength) { if (vars[ministep] == 2) { theBeh <- behChange(theBeh, egos[ministep], difs[ministep]) zlist <- c(zlist, list(theBeh)) } avz[ministep+1] <- mean(theBeh) } # What did we make: cbind(egos, vars, difs, avz) plot(avz) # Instead of avz, any function can be used. # zchain is the chain of z vectors, # but recorded only after ministeps for behavior: zchain <- matrix(unlist(zlist), nrow=length(startBeh), ncol=(1+sum(vars==2)), byrow=FALSE) # To see what was done, look at the beginning of the chain: vars[1:6] egos[1:6] difs[1:6] # This shows that in the first 6 ministeps, of which the first one # is a dummy, there were 2 behavioral ministeps, # of which the first one had a change of 0, # and the second one a change of +1 by actor 35. # This is reflected by row 35 of: zchain[,1:3] ############################################################################### #### Check by comparing the last state with the simulated behaviour #### #### as returned by siena07 #### ############################################################################### # The current value of theBeh is the last vector in zchain: all.equal(theBeh , zchain[,(1+sum(vars==2))]) # This is equal to the simulated value at the end of the period, # which is returned if returnDeps=TRUE. # Look at the help page for siena07 ?siena07 # to see what is returned if returnDeps=TRUE. # To check, simulate with the same algorithm - with the specified seed: sim2.ans <- siena07(sim.algorithm, data=mydata, effects=myeff, returnDeps=TRUE) # And indeed: all.equal(sim2.ans$sims[[runNumber]][[1]][[depVariable]][[period]] , theBeh)