############################################################################### #### WorkOnChains.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. # Written by Tom Snijders, with further help from Robert Krause, # version 13-06-2020 sim.algorithm <- sienaAlgorithmCreate(projname = 'prname', 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. theChain <- sim.ans$chain # The chain has the structure # chain[[run]][[group]][[period]][[ministep]] # (I'm not sure this is correct.) # For a single group, the [[group]] level is omitted. # For the meaning of the 13 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. # The stored ministeps give as "ReciRate" the reciprocal rates (expected values), # but not the actual time delays. # Therefore we need some extra effort to include time delays # with the correct distribution. # For this model these are i.i.d. and independent of the rest of the chain. # That is why they are not included in the chain: the time delays are not relevant, # only the parameters with which they were drawn. # So we can just generate them using the reciprocal rates. runNumber <- 1 (chainLength <- length(theChain[[runNumber]][[1]][[1]])) # Define the starting matrix: with imputed data. startMat <- ?? ## [data for wave 1] # Define the toggle function: toggle <- function(mat,i,j) { if (i != j) {mat[i,j] <- 1-mat[i,j]} mat } # Define the transitivity index transitivity <- function(mat){ m <- mat diag(m) <- 0 twop <- (m %*% m) # number of twopaths on each directed pair of nodes. diag(twop) <- 0 sum(m*twop)/sum(twop) # fraction of twopaths linked by ties. } # Note: this yields the same (as long as 0/0 is avoided) # as function gtrans in sna. # To calculate the intermediate networks # 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 first run, although 10 were generated. # Inspect the first three steps: theSteps[,1:3] # For the meaning of the 13 rows, see above. # Extract the relevant rows: egos <- as.numeric(theSteps[4,])+1 alters <- as.numeric(theSteps[5,])+1 egos <- c(0, egos) alters <- c(0, alters) invRates <- as.numeric(theSteps[7,]) # Prepare for calculations: avdegrees <- rep(NA, chainLength+1) trans.coefs <- rep(NA, chainLength+1) Delta <- rep(NA, chainLength+1) times <- rep(NA, chainLength+1) theMat <- startMat # Start calculations: times[1] <- 0.0 Delta[1] <- 0 trans.coefs[1] <- transitivity(theMat) avdegrees[1] <- sum(theMat)/dim(theMat)[1] # Pursue calculations: for (ministep in 1:chainLength) { theMat <- toggle(theMat, egos[ministep], alters[ministep]) if (egos[ministep] == alters[ministep]) { Delta[ministep] <- 0 } else { Delta[ministep] <- ifelse((theMat[egos[ministep], alters[ministep]] == 1), 1, -1) } avdegrees[ministep+1] <- sum(theMat)/dim(theMat)[1] trans.coefs[ministep+1] <- transitivity(theMat) # of course other functions can be added; or the entire theMat can be saved in a list times[ministep+1] <- times[ministep] + rexp(1, 1/invRates[ministep]) } # What did we make: cbind(times, egos, alters, Delta, avdegrees, trans.coefs) # And here is a function # to impute the next wave for dependent variable netName # the default number3=10 is what is perhaps the minimum number possible for siena07. imputeNetworkMissings <- function(seed, dat, eff, number3=10, netName){ toggle <- function(mat,i,j) { if (i != j) {mat[i,j] <- 1 - mat[i,j]} mat } alg <- sienaAlgorithmCreate(projname = 'imputeMissings', maxlike=TRUE, mult=20, seed=seed, nsub=0, n3=number3, simOnly=TRUE) if (dim(dat$depvars[[1]])[3] != 2){stop('More than 2 waves.')} wave1 <- dat$depvars[[1]][,,1] wave2 <- dat$depvars[[1]][,,2] diag(wave1) <- 0 diag(wave2) <- 0 if (any(is.na(wave1))){stop('Missings in wave1')} if (!any(is.na(wave2))){stop('No missings in wave2')} ans <- siena07(alg, data=dat, effects=eff, returnChains=TRUE) theChain <- ans$chain[[number3]][[1]][[1]] # take the last one chainLength <- length(theChain) theSteps <- matrix(unlist(theChain), 11, chainLength) # the following restriction is necessary # for cases where we have several dependent variables theSteps <- theSteps[,theSteps[3,]==netName] if (dim(theSteps)[2] == 0) { save(ans, file='saved_ans.RData') stop('Probably you used an incorrect netName: there are no steps with this name! But ans was saved.') } egos <- as.numeric(theSteps[4,])+1 alters <- as.numeric(theSteps[5,])+1 theMat <- wave1 for (ministep in 1:length(egos)) { theMat <- toggle(theMat, egos[ministep], alters[ministep]) } iseq <- (theMat == wave2) # is equal diag(iseq) <- TRUE iseq[is.na(iseq)] <- FALSE nisna2 <- (!is.na(wave2)) diag(nisna2) <- TRUE if (!isTRUE(all.equal(nisna2, iseq, check.attributes=FALSE))) { cat(all.equal(nisna2, iseq, check.attributes=FALSE),'\n') stop('Imputed matrix not for all observed values identical to wave 2.') } theMat }