Replicates and Variability

Steps followed:

• I produced 2000 replicates using ISA v6 for two scenarios: only influence and only selection.
• I computed average autocorrelation for each replicate (4 waves or measurements)
• I defined samples with no replacement of size n.
• I computed the coefficient of variation (CV = sd/mean) for each sample, then I plotted CV’s for each sample size.

rm(list=ls(all=TRUE))
#+ libraries
library(data.table)
library(RSienaTest)
library(igraph)
library(sna)
library(sdazar)
library(ggplot2)
gpath <- "/Users/sdaza/Desktop/model/"
folder_output <- "replicates"
source(paste0(gpath,"R/functions/auxiliary-functions.R"))
# compute network statistics
#+ set values
agents <- 200 # define number of agents to create network
iteration <- 1:2 # only two, influence and selection
lab_iteration <- c("Only Influence", "Only Selection")
nwaves <- c(6:9) # very important
max_rep <- 2000
#+ define paths
fpath <- paste0(gpath, "output/", folder_output) # files
fpath

[1] "/Users/sdaza/Desktop/model/output/replicates"

# get files
behavior.files <- list.files(path = fpath, pattern = "behavior")
b <- lapply(paste0(fpath, "/", behavior.files), fread, sep=",")
b <- rbindlist(b)
network.files <- list.files(path = fpath, pattern = "network")
n <- lapply(paste0(fpath, "/", network.files), fread, sep=",")
n <- rbindlist(n)
# categorize behavior variable
b$mybehavior <- cut(b$behavior, breaks = 10, labels = 1:10)
table(b$mybehavior, useNA = "ifany")

      1       2       3       4       5       6       7       8       9      10 
 560466  426152  580220  903571 1282580 1303592  939798  601516  436731  571574 

# loop through iterations
for (j in iteration) {
#+ create variables
lmoranb <- list()
rep <- 1:max_rep
#+ loop through replications
for (h in rep) {
# print(paste0(". . . iteration ", j, " replication ", h))
# extract data
tn <- n[iteration == j & replication == h]
tb <- b[iteration == j & replication == h]
# td <- d[iteration == j & replication == h]
#+ loop to create networks and get behavior
gnet <- list()
net <- list()
attr <- list()
# loop to get objects
for (i in 1:length(nwaves)) {
    t <- as.data.frame(tn[measurement == nwaves[i], 1:2])
    g <- graph.data.frame(t, vertices = 1:agents, directed = TRUE)
    m <- get.adjacency(g, sparse = FALSE)
    gnet[[i]] <- g
    net[[i]] <- m
    attr[[i]] <- tb[measurement == nwaves[i], list(id, mybehavior, radius, tendency)]
}
moranb <- NA
for (i in 1:length(nwaves)) {
  bh <- as.numeric(attr[[i]][, mybehavior])
  rd <- as.numeric(attr[[i]][, radius])
  moranb[i] <- Moran(bh, net[[i]])
}
lmoranb[[h]] <- mean(moranb)
} # end loop through replications
# coefficient of variation per replication
mb <- unlist(lmoranb)
CV <- function(x) (sd(x)/mean(x))
cvs <- list()
for (i in 2:500) { # size of samples
  s <- sample(rep, i, replace = FALSE)
  cvs[[i - 1]] <-  CV(mb[s])
}
cvs <- data.frame(cvs = unlist(cvs), n = 2:500)
print(ggplot(cvs, aes(y = cvs, x = n)) + 
  geom_point() + 
  labs(title = lab_iteration[j], x = "\nSample sizes", y =  "CV\n") + 
  geom_vline(xintercept = 30, linetype = 3))
} # end loop through iterations

Over 200 replicates there are not big gains in the reduction of CV’s behavior autocorrelation.