Appendix A — Permissions and Code for Chapter 2

This chapter contains the code used for the analysis in chapter 2 and the ODD (Overview, Design concepts, Details) protocol for the agent-based model ArchMatNet. The code is organized into three sections: R, Functions, and R Environment. The R section contains the code for the analysis, the Functions section contains the functions used in the analysis, and the R Environment section contains the R environment used to run the analysis. The final section contains the ODD protocol which is an enhanced description of the agent-based model.

R Code


library(tidyverse)
library(here)
library(rio)
library(ggpubr)
library(igraph)
library(ggnetwork)
library(tictoc)
# library(patchwork)
library(magick)
library(BayesFactor)
i_am("Analysis/Chap2Analysis.R")
source("Analysis/functions.R")

vars <- c(
 "uniqueTraits",
 "nTicks",
 "huntingGroupSize",
 "pHunting",
 "genderedMigration",
 "pNewAlliance",
 "maximumAlliesN",
 "pTrading",
 "visibility",
 "pLoss",
 "pLuckyLeap",
 "pVisiting",
 "learningMethod",
 "gradualEnvChange",
 "pEnvChange",
 "aggregationDuration",
 "transmissionRate",
 "pMigration",
 "nOfBands",
 "campPopulation",
 "potteryPrestige",
 "pNewVariant",
 "pNewObject",
 "aggregationFreq"
)

#### load results ####

library(arrow)
runs <- read_feather(here("data/modelOutput.parquet"))

#### look at range of experiments ####
expDF <- runs %>%
 select(expID) %>%
 group_by(expID) %>%
 count()

#### examine network comparisons ####

table_df <- tibble()

plotdf <- runs %>%
 filter(uniqueTraits == 1) %>%
 # filter(pTrading > 0) %>%
 # filter(nOfBands > 3) %>%
 select(contains("_corZscore")) %>%
 pivot_longer(contains("_corZscore"), names_to = "Network Comparison", values_to = "Correlation Comparison") %>%
 mutate(`Network Comparison` = gsub("_corZscore", "", `Network Comparison`)) %>%
 filter(!str_detect(`Network Comparison`, "Hunting")) %>%
 drop_na()

meanAll <- plotdf %>%
 group_by(`Network Comparison`) %>%
 summarize(mean = mean(`Correlation Comparison`, na.rm = T)) %>%
 arrange(mean)

table_df <- meanAll %>%
 mutate(metric = "Graph Correlation")

meanPotPoint <- plotdf %>%
 filter(`Network Comparison` != "Point_Pot") %>%
 filter(str_detect(`Network Comparison`, "Point|Pot")) %>%
 mutate(`Network Comparison` = case_when(`Network Comparison` %>% str_detect("Point") ~ "Point", TRUE ~ "Pot")) %>%
 group_by(`Network Comparison`) %>%
 summarize(mean = mean(`Correlation Comparison`)) %>%
 arrange(`Network Comparison`)

plotSum <- plotdf %>%
 group_by(`Network Comparison`) %>%
 summarize(mean = mean(`Correlation Comparison`, na.rm = T)) %>%
 arrange(desc(mean))
p <- plotdf %>%
 mutate(`Network Comparison` = fct_relevel(`Network Comparison`, plotSum$`Network Comparison`)) %>%
 ggplot(aes(x = `Network Comparison`, y = `Correlation Comparison`, fill = `Network Comparison`)) +
 geom_boxplot(alpha = .25, notch = T) +
 guides(fill = "none") +
 theme_bw() +
 xlab("") +
 theme(axis.text.x = element_text(angle = 45, hjust = 1))
p
# ggsave(here('figures/corComparisonZScoreBoxplot.png'),p,width = 7,height = 4.5, dpi = 600)
# ggsave(here('figures/corComparisonZScoreBoxplot.tiff'),p,width = 7,height = 4.5, compression = "lzw", dpi = 900)

plotdf <- runs %>%
 filter(uniqueTraits == 1) %>%
 select(contains("_QAP")) %>%
 pivot_longer(contains("_QAP"), names_to = "Network Comparison", values_to = "QAP Test") %>%
 mutate(`Network Comparison` = gsub("_QAP", "", `Network Comparison`)) %>%
 filter(!str_detect(`Network Comparison`, "Hunting")) %>%
 drop_na()

meanAll <- plotdf %>%
 group_by(`Network Comparison`) %>%
 summarize(mean = mean(`QAP Test`)) %>%
 arrange(mean)

table_df <- bind_rows(table_df, meanAll %>%
 mutate(metric = "QAP Test"))


meanPotPoint <- plotdf %>%
 filter(`Network Comparison` != "Point_Pot") %>%
 filter(str_detect(`Network Comparison`, "Point|Pot")) %>%
 mutate(`Network Comparison` = case_when(`Network Comparison` %>% str_detect("Point") ~ "Point", TRUE ~ "Pot")) %>%
 group_by(`Network Comparison`) %>%
 summarize(mean = mean(`QAP Test`)) %>%
 arrange(`Network Comparison`)
plotSum <- plotdf %>%
 group_by(`Network Comparison`) %>%
 summarize(mean = mean(`QAP Test`, na.rm = T)) %>%
 arrange(mean)
p <- plotdf %>%
 mutate(`Network Comparison` = fct_relevel(`Network Comparison`, plotSum$`Network Comparison`)) %>%
 ggplot(aes(x = `Network Comparison`, y = `QAP Test`)) +
 geom_point(position = position_jitter(), alpha = .1) +
 theme_bw() +
 xlab("") +
 theme(axis.text.x = element_text(angle = 45, hjust = 1))
p

# ggsave(here('figures/QAPTest.tiff'),p,width = 7,height = 4.5, compression = "lzw", dpi = 900)
# browseURL(here('figures/QAPTest.tiff'))

# table_df  %>%
# pivot_wider(names_from = metric, values_from = mean) %>%
# write_csv("tables/MetricResults.csv")

#### explore variables ####

experiments <- runs$expID %>%
 unique() %>%
 sort()
experiments

# unique traits

i <- 24
variables_df <- runs %>%
 select(expID, any_of(vars), any_of(contains("_corZscore"))) %>%
 pivot_longer(any_of(contains("_corZscore"))) %>%
 filter(!str_detect(name, "Hunting")) %>%
 filter(str_detect(name, "Pot|Point")) %>%
 pivot_longer(any_of(vars), names_to = "Variable", values_to = "parameter") %>%
 mutate(parameter = factor(parameter)) %>%
 group_by(name, Variable, parameter) %>%
 summarize(cor = mean(value, na.rm = T), count = n())

visiting_point_df <- variables_df %>%
 filter(name == "Visiting_Point_corZscore") %>%
 filter(Variable == vars[i])

visiting_point_df

runs %>%
 select(uniqueTraits, nTicks, Visiting_Point_corZscore) %>%
 mutate_at(vars(-Visiting_Point_corZscore), factor) %>%
 ggplot(aes(x = Visiting_Point_corZscore, group = nTicks, color = nTicks)) +
 geom_boxplot() +
 theme_minimal() +
 facet_wrap(~uniqueTraits, scales = "free_y") +
 ggtitle("uniqueTraits")

variables_df <- runs %>%
 filter(uniqueTraits == 1) %>%
 select(-uniqueTraits) %>%
 select(expID, any_of(vars), any_of(contains("_corZscore"))) %>%
 pivot_longer(any_of(contains("_corZscore"))) %>%
 filter(!str_detect(name, "Hunting")) %>%
 filter(str_detect(name, "Pot|Point")) %>%
 filter(name == "Visiting_Pot_corZscore") %>%
 janitor::remove_constant()

library(mgcv)

dvars <- vars[which(vars %in% names(variables_df))]

fmla <- as.formula("value ~ nTicks + huntingGroupSize + pHunting +
 factor(genderedMigration) + pNewAlliance + maximumAlliesN +
 pTrading + factor(visibility) + pLoss + pLuckyLeap +
 pVisiting + factor(learningMethod) + factor(gradualEnvChange) +
 pEnvChange + aggregationDuration + transmissionRate +
 pMigration + nOfBands + campPopulation +
 factor(potteryPrestige) + pNewVariant + pNewObject +
 aggregationFreq")
gam_mdl <- gam(fmla, data = variables_df)
summary(gam_mdl)
gam_summary <- summary(gam_mdl)

# Extract parametric and smooth terms
parametric_df <- as.data.frame(gam_summary$p.table)
smooth_df <- as.data.frame(gam_summary$s.table)

# Add type column and combine tables
parametric_df$type <- "Parametric"

# Format table
names(parametric_df) <- c("Estimate", "Std. Error", "t-value", "p-value", "Type")

parametric_df <- rownames_to_column(parametric_df, "Parameter")

parametric_df <- parametric_df %>%
 mutate_if(is.numeric, round, 2)

# write_csv(parametric_df,here('tables/ParametricEstimates.csv'))

# variables that don't make a difference
xvars <- c("pNewAlliance", "maximumAlliesN", "pLoss", "pEnvChange", "aggregationDuration", "potteryPrestige", "pNewVariant", "transmissionRate", "pMigration", "nOfBands", "pNewVariant", "pLuckyLeap", "pHunting")

dvars <- vars[which(!vars %in% xvars)]
dvars <- dvars[which(dvars %in% names(variables_df))]

# fmla = as.formula(paste("value ~ ",paste(dvars,collapse = "+")))
fmla <- as.formula("value ~ nTicks +
 huntingGroupSize +
 factor(genderedMigration) +
 pTrading + factor(visibility) +
 pVisiting + factor(learningMethod) + factor(gradualEnvChange) + campPopulation + pNewObject +
 aggregationFreq")

gam_mdl2 <- gam(fmla, data = variables_df)
summary(gam_mdl2)

xvars2 <- c("pNewAlliance", "maximumAlliesN", "pLoss", "pEnvChange", "aggregationDuration", "potteryPrestige", "pNewVariant", "transmissionRate", "pMigration", "nOfBands", "pNewVariant", "pLuckyLeap", "genderedMigration", "pHunting", "pTrading", "gradualEnvChange", "huntingGroupSize", "pNewObject", "pVisiting", "campPopulation")

dvars <- vars[which(!vars %in% xvars2)]
dvars <- dvars[which(dvars %in% names(variables_df))]

# fmla = as.formula(paste("value ~ ",paste("factor(",dvars,")",collapse = "+")))
fmla <- as.formula("value ~
 nTicks +
 factor(visibility) +
 factor(learningMethod) +
 aggregationFreq")

gam_mdl3 <- gam(fmla, data = variables_df)
summary(gam_mdl3)

gam_nTicks <- gam(value ~ factor(nTicks), data = variables_df)
gam_visibility <- gam(value ~ factor(visibility), data = variables_df)
gam_pVisiting <- gam(value ~ pVisiting, data = variables_df)
gam_learningMethod <- gam(value ~ factor(learningMethod), data = variables_df)
gam_aggregationFreq <- gam(value ~ aggregationFreq, data = variables_df)
gam_campPopulation <- gam(value ~ campPopulation, data = variables_df)
summary(gam_nTicks)
summary(gam_visibility)
summary(gam_pVisiting)
summary(gam_learningMethod)
summary(gam_aggregationFreq)
summary(gam_campPopulation)


variables_df %>%
 select(-any_of(xvars2)) %>%
 pivot_longer(any_of(vars), names_to = "Variable", values_to = "parameter") %>%
 mutate(parameter = factor(parameter)) %>%
 rename(cor = value) %>%
 ggplot(aes(x = parameter, y = cor, color = Variable)) +
 geom_boxplot() +
 guides(color = "none") +
 theme_minimal() +
 facet_wrap(~Variable, scales = "free_x")

# ggsave(here('figures/VariableCorrelations.tiff'),width = 7,height = 4.5, compression = "lzw", dpi = 900)
# ggsave(here('figures/VariableCorrelations.png'),width = 7,height = 4.5, dpi = 150)
# browseURL(here('figures/VariableCorrelations.tiff'))

#### look at individual networks ####
tictoc::tic()
g <- runs %>%
 filter(nOfBands == 5) %>%
 filter(uniqueTraits == 0) %>%
 filter(!is.na(wordMatrixToRowListCampVisitingAdjacency)) %>%
 slice_sample(n = 100) %>%
 mutate_all(as.character) %>%
 select(-runNumber) %>%
 rowid_to_column("runNumber") %>%
 pivot_longer(-runNumber,
   names_to = "property",
   values_to = "value"
 ) %>%
 filter(property %in% vars | str_detect(property, "wordMatrixToRowListCamp|runNumber")) %>%
 mutate(property = str_remove_all(
   property,
   "wordMatrixToRowListCamp|Adjacency"
 )) %>%
 pivot_wider(names_from = property, values_from = value) %>%
 pivot_longer(Visiting:Pot, names_to = "property", values_to = "adj") %>%
 filter(!is.na(adj)) %>%
 mutate(property = tolower(property)) %>%
 # rowwise() %>%
 mutate(adj = formatNLAdjTables(adj)) %>%
 pivot_wider(names_from = property, values_from = adj)
# ungroup()
tictoc::toc()
bands <- g$point[1][[1]] %>%
 select(name = from) %>%
 distinct_all() %>%
 bind_rows(g$point[1][[1]] %>%
   select(name = to)) %>%
 distinct_all() %>%
 mutate(band = sort(rep(0:4, 3)))

getGraphDF <- function(network, i = 1) {
 cutoff <- quantile(g[[network]][[i]]$weight, .5)
 gdf <- g[[network]][[i]] %>%
   filter(weight > cutoff)
 n <- graph_from_data_frame(gdf, directed = F)
 # summary(n)
 plotdf <- ggnetwork(n) %>%
   select(-band) %>%
   left_join(bands) %>%
   mutate(
     band = factor(band),
     network = !!network
   )


 return(plotdf)
}

variables <- c("visiting", "learning", "trading", "point", "pot")
combined_df <- map_df(variables, getGraphDF) %>%
 mutate(network = factor(network, levels = variables))

p <- ggplot(combined_df, aes(x = x, y = y, xend = xend, yend = yend)) +
 geom_edges(aes(color = weight, alpha = weight), linewidth = 1) +
 geom_nodelabel(aes(fill = band, label = name), size = 3, color = "white") +
 scale_fill_viridis_d() +
 theme_void() +
 coord_equal(ratio = 1 / 2, clip = "off", expand = T) +
 theme(
   legend.position = "none",
   plot.background = element_rect(colour = "#404040", fill = NA, linewidth = .5),
   plot.title = element_text(hjust = 0.5),
   plot.margin = unit(c(2, 2, 2, 2), "mm"),
   strip.placement = "outside",
   strip.text = element_text(margin = margin(b = 10))
 ) +
 facet_wrap(~network)

p

# ggsave("MaterialCultureNetworksPaper\\figures\\SameTraits.png",dpi = 600, width = 8, height = 4.5, bg = "white")
# fn = "MaterialCultureNetworksPaper\\figures\\SameTraits.tiff"
# ggsave(fn,dpi = 900, width = 7, height = 5, bg = "white", compression = "lzw")

# img = image_read(fn)
# img = image_trim(img)
# img = image_write(img,fn)
#
# image1 = image_read("MaterialCultureNetworksPaper\\figures\\SameTraits.tiff")
# image2 = image_read("MaterialCultureNetworksPaper\\figures\\UniqueTraits.tiff")

# image1_annotated <- image_annotate(image1, text = "a",
#                                    gravity = "northwest",
#                                    location = "+20+20",
#                                    size = 200, color = "black")
#
# image2_annotated <- image_annotate(image2, text = "b",
#                                    gravity = "northwest",
#                                    location = "+20+20",
#                                    size = 200, color = "black")

# stacked_images <- image_append(c(image1_annotated, image2_annotated), stack = TRUE)
# image_write(stacked_images, "MaterialCultureNetworksPaper\\figures\\combined_networks.tiff")

# g1 = getGraph("visiting",2)
# g2 = getGraph("learning",2)
# g3 = getGraph("trading",2)
# g4 = getGraph("point",2)
# g5 = getGraph("pot",2)
#
# g1 + g2 + g3 + g4 + g5
# ggsave("MaterialCultureNetworksPaper\\figures\\SameTraits2.png",dpi = 600, width = 8, height = 4.5, bg = "white")

g <- runs %>%
 filter(nOfBands > 3) %>%
 filter(uniqueTraits == 1) %>%
 filter(!is.na(wordMatrixToRowListCampVisitingAdjacency)) %>%
 slice_sample(n = 100) %>%
 mutate_all(as.character) %>%
 select(-runNumber) %>%
 rowid_to_column("runNumber") %>%
 pivot_longer(-runNumber,
   names_to = "property",
   values_to = "value"
 ) %>%
 filter(property %in% vars | str_detect(property, "wordMatrixToRowListCamp|runNumber")) %>%
 mutate(property = str_remove_all(
   property,
   "wordMatrixToRowListCamp|Adjacency"
 )) %>%
 pivot_wider(names_from = property, values_from = value) %>%
 pivot_longer(Visiting:Pot, names_to = "property", values_to = "adj") %>%
 filter(!is.na(adj)) %>%
 mutate(property = tolower(property)) %>%
 # rowwise() %>%
 mutate(adj = formatNLAdjTables(adj)) %>%
 pivot_wider(names_from = property, values_from = adj)

variables <- c("visiting", "learning", "trading", "point", "pot")
combined_df <- map_df(variables, getGraphDF) %>%
 mutate(network = factor(network, levels = variables))

p <- ggplot(combined_df, aes(x = x, y = y, xend = xend, yend = yend)) +
 geom_edges(aes(color = weight, alpha = weight), linewidth = 1) +
 geom_nodelabel(aes(fill = band, label = name), size = 3, color = "white") +
 scale_fill_viridis_d() +
 theme_void() +
 coord_equal(ratio = 1 / 2, clip = "off", expand = T) +
 theme(
   legend.position = "none",
   plot.background = element_rect(colour = "#404040", fill = NA, linewidth = .5),
   plot.title = element_text(hjust = 0.5),
   plot.margin = unit(c(2, 2, 2, 2), "mm"),
   strip.placement = "outside",
   strip.text = element_text(margin = margin(b = 10))
 ) +
 facet_wrap(~network)

p

# fn = "MaterialCultureNetworksPaper\\figures\\UniqueTraits.tiff"
# ggsave(fn,dpi = 900, width = 7, height = 5, bg = "white", compression = "lzw")
#
# img = image_read(fn)
# img = image_trim(img)
# img = image_write(img,fn)

#### centrality ####

largeNetworks <- runs %>%
 filter(nOfBands > 3) %>%
 filter(uniqueTraits == 1)
results_point_e <- list()
results_pot_e <- list()
for (i in 1:nrow(largeNetworks)) {
 print(paste(i, "of", nrow(largeNetworks)))
 x <- largeNetworks[i, ]

 n1 <- matrix2zscore(x[["wordMatrixToRowListCampVisitingAdjacency"]])
 diag(n1) <- 0
 ec1 <- sna::evcent(n1)
 n2 <- matrix2zscore(x[["wordMatrixToRowListCampPointAdjacency"]])
 diag(n2) <- 0
 ec2 <- sna::evcent(n2)
 n3 <- matrix2zscore(x[["wordMatrixToRowListCampPotAdjacency"]])
 diag(n3) <- 0
 ec3 <- sna::evcent(n3)
 r1 <- cor(ec2, ec1) %>% round(., 2)
 r2 <- cor(ec3, ec1) %>% round(., 2)
 results_point_e[[i]] <- r1
 results_pot_e[[i]] <- r2
}


results_point_e <- unlist(results_point_e)
results_pot_e <- unlist(results_pot_e)
hist(results_point_e)
hist(results_pot_e)
results_point_e %>% summary()
results_pot_e %>% summary()

cor(largeNetworks$nOfBands, results_point_e)

plot(largeNetworks$nOfBands, results_point_e)

results_point_d <- list()
results_pot_d <- list()
for (i in 1:nrow(largeNetworks)) {
 print(paste(i, "of", nrow(largeNetworks)))
 x <- largeNetworks[i, ]

 n1 <- matrix2zscore(x[["wordMatrixToRowListCampVisitingAdjacency"]])
 diag(n1) <- 0
 dc1 <- sna::degree(n1)
 n2 <- matrix2zscore(x[["wordMatrixToRowListCampPointAdjacency"]])
 diag(n2) <- 0
 dc2 <- sna::degree(n2)
 n3 <- matrix2zscore(x[["wordMatrixToRowListCampPotAdjacency"]])
 diag(n3) <- 0
 dc3 <- sna::degree(n3)
 r1 <- cor(dc2, dc1) %>% round(., 2)
 r2 <- cor(dc3, dc1) %>% round(., 2)
 results_point_d[[i]] <- r1
 results_pot_d[[i]] <- r2
}

results_point_d <- unlist(results_point_d)
results_pot_d <- unlist(results_pot_d)
hist(results_point_d)
hist(results_pot_d)
results_point_d %>% summary()
results_pot_d %>% summary()

results_point_b <- list()
results_pot_b <- list()
for (i in 1:nrow(largeNetworks)) {
 print(paste(i, "of", nrow(largeNetworks)))
 x <- largeNetworks[i, ]

 n1 <- matrix2zscore(x[["wordMatrixToRowListCampVisitingAdjacency"]])
 diag(n1) <- 0
 bc1 <- sna::betweenness(n1)
 n2 <- matrix2zscore(x[["wordMatrixToRowListCampPointAdjacency"]])
 diag(n2) <- 0
 bc2 <- sna::betweenness(n2)
 n3 <- matrix2zscore(x[["wordMatrixToRowListCampPotAdjacency"]])
 diag(n3) <- 0
 bc3 <- sna::betweenness(n3)
 r1 <- cor(bc2, bc1) %>% round(., 2)
 r2 <- cor(bc3, bc1) %>% round(., 2)
 results_point_b[[i]] <- r1
 results_pot_b[[i]] <- r2
}

results_point_b <- unlist(results_point_b)
results_pot_b <- unlist(results_pot_b)
hist(results_point_b)
hist(results_pot_b)
results_point_b %>% summary()
results_pot_b %>% summary()

results_point_c <- list()
results_pot_c <- list()
for (i in 1:nrow(largeNetworks)) {
 print(paste(i, "of", nrow(largeNetworks)))
 x <- largeNetworks[i, ]

 n1 <- matrix2zscore(x[["wordMatrixToRowListCampVisitingAdjacency"]])
 diag(n1) <- 0
 bc1 <- sna::closeness(n1)
 n2 <- matrix2zscore(x[["wordMatrixToRowListCampPointAdjacency"]])
 diag(n2) <- 0
 bc2 <- sna::closeness(n2)
 n3 <- matrix2zscore(x[["wordMatrixToRowListCampPotAdjacency"]])
 diag(n3) <- 0
 bc3 <- sna::closeness(n3)
 r1 <- cor(bc2, bc1) %>% round(., 2)
 r2 <- cor(bc3, bc1) %>% round(., 2)
 results_point_c[[i]] <- r1
 results_pot_c[[i]] <- r2

 # ng = n1
 # ng[ng > .75] = 1
 # ng[ng <= .75] = 0
 # g = ng %>% graph_from_adjacency_matrix(mode = "undirected", weighted = F)
 # degree_distribution(g)
}

results_point_c <- unlist(results_point_c)
results_pot_c <- unlist(results_pot_c)
hist(results_point_c)
hist(results_pot_c)
results_point_c %>% summary()
results_pot_c %>% summary()

# examine bad network correlation
library(mgcv)
points_df <- runs %>%
 filter(uniqueTraits == 1) %>%
 select(expID, any_of(vars), any_of(contains("_corZscore"))) %>%
 pivot_longer(any_of(contains("_corZscore"))) %>%
 filter(name == "Visiting_Point_corZscore") %>%
 pivot_longer(vars, names_to = "variable", values_to = "parameter") %>%
 group_by(variable, parameter) %>%
 summarize(cor = mean(value, na.rm = T)) %>%
 mutate(parameter = as.character(parameter))

run_sample <- runs %>%
 filter(uniqueTraits == 1) %>%
 filter(pTrading > 0) %>%
 filter(learningMethod == 50) %>%
 filter(Visiting_Point_corZscore < 0) %>%
 slice_sample(n = 1)
run_sample$experiment
run_sample$runNumber
run_sample$nTicks

default_vars <- defaults %>%
 mutate(default = TRUE) %>%
 select(parameters = Parameter, default_values = Default, default)


#### ERGM ####

library(statnet)
library(Rglpk)
library(ergm)
library(dplyr)
library(igraph)
# compare visiting and other networks to pot and point
# compare nonunique traits to unique traits
# statistics to include mcmc.diagnostics - check simulated distributions have a peak at the estimate and goodness of fit diagnostics
# are the same generative processes responsible for all networks? If we just had the material culture then would we still come to the same conclusions?

set.seed(1010)

run_ergm <- function(i, network, formula_type) {
 formula <- switch(formula_type,
   "edges" = ns ~ edges + nodematch("band"),
   "triangle" = ns ~ edges + triangle + nodematch("band"),
   "threetrail" = ns ~ edges + triangle + threetrail + nodematch("band"),
   "gwesp" = ns ~ edges + gwesp(0.25, fixed = TRUE) + threetrail + nodematch("band"),
   stop("Unknown formula type")
 )

 library(statnet)
 library(igraph)
 library(dplyr)
 g <- readRDS("data/graphSample.Rds")
 tryCatch(
   {
     cutoff <- quantile(g[[network]][[i]]$weight, .75)
     gdf <- g[[network]][[i]] %>%
       filter(weight > cutoff)
     n <- graph_from_data_frame(gdf, directed = F)

     nodeList <- gdf %>%
       select(id = from, band) %>%
       distinct_all()
     ns <- intergraph::asNetwork(n)
     ns %v% "band" <- nodeList$band
     controlVals <- control.ergm(
       parallel = 14,
       seed = 1010,
       MCMLE.confidence = 0.9,
       MPLE.maxit = 1000
     )
     x <- tryCatch(
       ergm(formula, control = controlVals),
       error = function(e) {
         print(e)
         return(NA)
       }
     )
     e <- tryCatch(
       summary(x),
       error = function(e) {
         return(NA)
       }
     )
     gf <- tryCatch(
       gof(x),
       error = function(e) {
         return(NA)
       }
     )
     f <- tryCatch(
       x$formula %>% as.character() %>% paste(collapse = " "),
       error = function(e) {
         return(NA)
       }
     )
     aic <- tryCatch(
       e$aic,
       error = function(e) {
         return(NA)
       }
     )
     gf_pval <- tryCatch(
       gf$pval.model[1, 5],
       error = function(e) {
         return(NA)
       }
     )
     edges <- tryCatch(
       e$coefficients["edges", "Pr(>|z|)"],
       error = function(e) {
         return(NA)
       }
     )
     result <- tibble(
       run = i,
       network = network,
       formula = f,
       aic = aic,
       gf_pval = gf_pval,
       edges = edges
     )
     saveRDS(
       result,
       paste0(
         "data/ergm_results_",
         i,
         "_",
         network,
         "_",
         formula_type,
         ".rds"
       )
     )
     return(result)
   },
   error = function(e) {
     return(NA)
   }
 )
}

run_and_monitor_processes <- function(bg, timeout = 120) {
 start_times <- Sys.time()

 # Monitor processes
 while (TRUE) {
   Sys.sleep(1) # Check every second

   # Check the status of each process
   for (name in names(bg)) {
     if (bg[[name]]$is_alive()) {
       run_time <- as.numeric(difftime(Sys.time(), start_times, units = "secs"))
       if (shiny::isTruthy(run_time > timeout)) {
         bg[[name]]$kill() # Kill process if it exceeds timeout
         cat("Process", name, "killed after", run_time, "seconds\n")
       }
     }
   }

   # Exit the loop if all processes are done
   if (all(!sapply(bg, function(p) p$is_alive()))) {
     break
   }
 }
}

# run_ergm(1, "trading", "edges")
library(callr)

bg <<- list()

for (i in 1:10) {
 print(i)

 bg <- list()

 for (network in c("trading", "learning", "visiting", "point", "pot")) {
   print(network)
   bg <- list()

   for (formula_type in c("edges", "triangle", "threetrail", "gwesp")) {
     l <- list(i = i, network = network, formula_type = formula_type)
     bg[[paste0(i, "-", network, "-", formula_type)]] <- r_bg(run_ergm, args = l, supervise = FALSE)
   }
   run_and_monitor_processes(bg, timeout = 120)
 }
}

ls <- list.files("data/", pattern = "ergm_results", full.names = T)
ls <- ls[which(!str_detect(ls, ".csv"))]

results <- lapply(ls, readRDS) %>% bind_rows()


results_summary <- results %>%
 group_by(network, formula) %>%
 summarize(mean_gf_pval = mean(gf_pval, na.rm = T), mean_edges = mean(edges, na.rm = T), mean_aic = mean(aic, na.rm = T), count = n()) %>%
 ungroup() %>%
 filter(count > 15)

results_table <- results_summary %>%
 select(network, formula, gof = mean_gf_pval, aic = mean_aic) %>%
 filter(network %in% c("visiting", "point")) %>%
 mutate(formula = case_when(
   str_detect(formula, "gwesp") ~ "edge + gwesp + threetrail + band",
   str_detect(formula, "threetrail") ~ "edge + triangle + threetrail + band",
   TRUE ~ "edge + band"
 )) %>%
 arrange(network, formula)
results_table
results_table %>% write_csv("tables/ergm_results_table.csv")

results_summary2 <- results %>%
 group_by(formula) %>%
 summarize(mean_gf_pval = mean(gf_pval), mean_edges = mean(edges), mean_aic = mean(aic))
results_summary2
# graphical representation of how an ergm works

library(igraph)
set.seed(1010)
g <- readRDS("data/graphSample.Rds")
i <- 10
network <- "point"
cutoff <- quantile(g[[network]][[i]]$weight, .75)
gdf <- g[[network]][[i]] %>%
 filter(weight > cutoff)
n <- graph_from_data_frame(gdf, directed = F)

# nodeList = g[[network]][[i]] %>% select(id = from, band) %>%
# distinct_all()
nodeList <- structure(
 list(id = c(
   "5", "6", "7", "8", "9", "10", "11", "12",
   "13", "14", "15", "16", "17", "18", "19"
 ), band = c(
   "0", "0", "0",
   "1", "1", "1", "2", "2", "2", "3", "3", "3", "4", "4", "4"
 )),
 row.names = c(
   NA,
   -15L
 ), class = c("tbl_df", "tbl", "data.frame")
)
g_df <- tibble(id = V(n)$name) %>%
 left_join(nodeList)
V(n)$band <- g_df$band
band_colors <- c("0" = "red", "1" = "blue", "2" = "green", "3" = "orange", "4" = "purple")
V(n)$color <- band_colors[V(n)$band]
# png("figures/ergm_example.png", width = 800, height = 800, bg = "black")
# plot(n, main = "", vertex.label.color = "black")
# title(main = "Observed Network", col.main = "white")
# dev.off()
# g_random <- rewire(n, keeping_degseq(niter = vcount(n) * 10))
# png("figures/ergm_example_random.png", width = 800, height = 800, bg = "black")
# plot(g_random, main = "", vertex.label.color = "black")
# title(main = "Random Network", col.main = "white")
# dev.off()
# g <- make_ring(15) %>%
#   add_edges(c(1,5, 2,6, 3,7, 4,8, 5,9, 6,10, 7,11, 8,12, 9,13, 10,14, 11,15))
# plot(g)
ns <- intergraph::asNetwork(n)
# plot(ns)
formula_type <- "gwesp"
formula <- switch(formula_type,
 "edges" = ns ~ edges + nodematch("band"),
 "triangle" = ns ~ edges + triangle + nodematch("band"),
 "threetrail" = ns ~ edges + triangle + threetrail + nodematch("band"),
 "gwesp" = ns ~ edges + gwesp(0.25, fixed = TRUE) + threetrail + nodematch("band"),
 stop("Unknown formula type")
)

random_ergm <- ergm(formula)
summary(random_ergm)
pgof <- gof(random_ergm)
pdf(glue::glue("figures/gof_{network}_{formula_type}.pdf"))
plot(pgof)
dev.off()
# test probability that random networks are correlated

g <- readRDS("data/graphSample.Rds")
i <- sample(1:100, 1)
i
network <- "visiting"
cutoff <- quantile(g[[network]][[i]]$weight, .75)
gdf <- g[[network]][[i]] %>%
 filter(weight > cutoff)
n <- graph_from_data_frame(gdf, directed = F)

nodeList <- gdf %>%
 select(id = from, band) %>%
 distinct_all()
ns <- intergraph::asNetwork(n)
ns %v% "band" <- nodeList$band
controlVals <- control.ergm(
 parallel = 14,
 seed = 1010,
 MCMLE.confidence = 0.9,
 MPLE.maxit = 1000
)

ergm_model <- ergm(ns ~ edges + triangle + threetrail + nodematch("band"), control = controlVals)
# ns_visiting = ns
nsim <- simulate(ergm_model, nsim = 10000)
gcor_values <- numeric(length = length(nsim))
# ns_matrix = as.matrix(ns)
ns_matrix <- as.matrix(ns_visiting)


# Loop through each simulated network and calculate gcor
for (i in seq_along(nsim)) {
 nsim_matrix <- as.matrix(nsim[[i]]) # Convert simulated network to matrix
 ns_row_names <- rownames(ns_matrix)
 ns_col_names <- colnames(ns_matrix)
 nsim_row_names <- rownames(nsim_matrix)
 nsim_col_names <- colnames(nsim_matrix)

 # Union of row and column names (all possible nodes in both matrices)
 all_row_names <- union(ns_row_names, nsim_row_names)
 all_col_names <- union(ns_col_names, nsim_col_names)

 # Initialize matrices to the larger set of row and column names, filled with zeros
 new_ns_matrix <- matrix(0,
   nrow = length(all_row_names), ncol = length(all_col_names),
   dimnames = list(all_row_names, all_col_names)
 )
 new_nsim_matrix <- matrix(0,
   nrow = length(all_row_names), ncol = length(all_col_names),
   dimnames = list(all_row_names, all_col_names)
 )

 # Fill the new matrices with the values from the original matrices
 new_ns_matrix[ns_row_names, ns_col_names] <- ns_matrix
 new_nsim_matrix[nsim_row_names, nsim_col_names] <- nsim_matrix
 gcor_values[i] <- gcor(new_ns_matrix, new_nsim_matrix) # Calculate gcor
}

hist(gcor_values)
gcor_values %>% summary()



network <- "visiting"
cutoff <- quantile(g[[network]][[i]]$weight, .75)
gdf <- g[[network]][[i]] %>%
 filter(weight > cutoff)
n <- graph_from_data_frame(gdf, directed = F)

nodeList <- gdf %>%
 select(id = from, band) %>%
 distinct_all()
ns <- intergraph::asNetwork(n)
ns %v% "band" <- nodeList$band
controlVals <- control.ergm(
 parallel = 14,
 seed = 1010,
 MCMLE.confidence = 0.9,
 MPLE.maxit = 1000
)

ergm_model <- ergm(ns ~ edges + triangle + threetrail + nodematch("band"), control = controlVals)
nsim <- simulate(ergm_model, nsim = 1000)
gcor_values <- numeric(length = length(nsim))
ns_matrix <- as.matrix(ns)

# Loop through each simulated network and calculate gcor
for (i in seq_along(nsim)) {
 nsim_matrix <- as.matrix(nsim[[i]]) # Convert simulated network to matrix
 gcor_values[i] <- gcor(ns_matrix, nsim_matrix) # Calculate gcor
}

# hist(gcor_values)
# gcor_values %>% summary()


#### export session info  ####
# Get the R version
r_version <- paste0("R version=", paste(R.version$major, R.version$minor, sep = "."))

# Get the loaded packages and their versions
loaded_packages <- sessionInfo()$otherPkgs
loaded_versions <- sapply(loaded_packages, function(x) paste(x$Package, x$Version, sep = "="))

# Write the R version and package versions to the file
writeLines(c(r_version, loaded_versions), con = "analysis/r_environment.qmd")

Functions

# custom functions
library(tictoc)
library(magrittr)
library(ggnetwork)
library(ggthemes)

#' Load all data by specifying one complete name or the partial date file name
#'
#' @param fileName
#'
#' @return
#' @export
#'
#' @examples
loadRun <- function(fileName) {
  if (file.exists(here(stringr::str_glue("results/{fileName}")))) {
    fileName %<>% str_remove_all("run_[A-z]*_")
  }
  map(
    list.files(here("results"),
      pattern = fileName,
      full.names = T
    ) %>% .[str_detect(., fileName)],
    ~
      read_csv(.x, show_col_types = F)
  ) %>% set_names(c(
    "assemblages",
    "links",
    "nodes",
    "variables"
  ))
}

# x = runResults$assemblages$pointAssemblage[[7]]

#' Convert NetLogo list format into R list
bracket2list <- function(x) {
  if (str_detect(x, "[0-9]")) {
    x %>%
      str_remove_all("\\]|\\[") %>%
      str_split_fixed(" ", n = Inf) %>%
      matrix(ncol = 2, byrow = T) %>%
      as_tibble(.name_repair = "minimal") %>%
      setNames(c("tick", "tool"))
  } else {
    NA
  }
}


#' Calculate manhattan similarity
#'
#' @param analysis dataframe formatted through formatAssemblage function
#' @param type pot or point
#' @param group band or site
#' @param cutoff quantile cutoff to keep links
#' @param similarity either 'BR' for Brainerd-Robinson or 'WJ' for weighted Jaccard
#'
#' @return dataframe of nodes and dataframe of links
#' @export
#'
#' @examples
calcDistanceTools <- function(analysis,
                              toolType = "pot",
                              group = "band",
                              cutoff = 0,
                              similarity = "WJ") {
  tmp <- analysis %>%
    filter(type == toolType)
  m <- table(tmp[[group]], tmp$tool) %>%
    prop.table(1) %>%
    {
      . * 100
    }

  # calculate similarity
  if (similarity == "BR") {
    g <- dist(m, "manhattan") %>%
      as.matrix()
    g <- 1 - (g / 200)
    diag(g) <- 0
    # g[] = round(g, 2)
  } else {
    g <- weightedJaccard(m)
    # g[] = round(g, 2)
  }

  links <- g %>%
    as.data.frame() %>%
    rownames_to_column("from") %>%
    pivot_longer(-from,
      names_to = "to",
      values_to = "weight"
    ) %>%
    # remove duplicate links
    dplyr::rowwise() %>%
    mutate(combo = paste(min(from, to), max(from, to))) %>%
    arrange(desc(to)) %>%
    distinct(combo, .keep_all = T) %>%
    select(-combo)
  links %<>%
    filter(weight > 0)
  cutoff <- quantile(links$weight, cutoff)
  links %<>%
    filter(weight >= cutoff) %>%
    mutate_at(vars(weight), round, 2)
  if (group != "band" & "band" %in% names(analysis)) {
    nodes <- links %>%
      select(from, to) %>%
      pivot_longer(c(from, to), values_to = "name", names_to = "var") %>%
      select(name) %>%
      arrange(name) %>%
      distinct_all() %>%
      full_join(
        analysis %>%
          select(name = site, band) %>%
          mutate_at(vars(name), as.character),
        by = "name"
      ) %>%
      mutate(camp = name) %>%
      distinct_all() %>%
      mutate_at(vars(band), factor)
  } else {
    nodes <- analysis %>%
      ungroup() %>%
      select(any_of(group)) %>%
      mutate(name = !!as.name(group)) %>%
      distinct_all()
  }

  return(list(links = links, nodes = nodes))
}

# function to remove labels from names
rmLabel <- function(x) {
  map_chr(
    x,
    ~ str_extract_all(.x, "[0-9]") %>%
      unlist() %>%
      paste(collapse = "")
  )
}

# links = runResults$links
# nodes = runResults$nodes
calcDistanceLinks <- function(links,
                              nodes,
                              colS = "all",
                              group,
                              cutoff = 0) {
  if (colS == "all") {
    colS <- c("learning", "hunting", "visiting", "trading")
  }
  df <- links %>%
    select(end1, end2, all_of(colS)) %>%
    mutate_at(vars(end1, end2), rmLabel) %>%
    left_join(
      nodes %>%
        select(end1 = who, band, camp) %>%
        mutate(end1 = end1 %>% as.character()),
      by = "end1"
    ) %>%
    left_join(
      nodes %>%
        select(end2 = who, band, camp) %>%
        mutate(end2 = end2 %>% as.character()),
      by = "end2",
      suffix = c("", "2")
    )
  linksMod <- df %>%
    select(group, any_of(str_glue("{group}2")), all_of(colS)) %>%
    # create unique paths
    mutate(path = map2_chr(
      !!as.name(group),
      !!as.name(str_glue("{group}2")),
      ~ sort(c(.x, .y)) %>% paste(collapse = "_")
    )) %>%
    group_by(path) %>%
    summarize_at(vars(all_of(colS)), sum) %>%
    ungroup() %>%
    separate(path, into = c(group, str_glue("{group}2")), sep = "_") %>%
    mutate_at(vars(group, str_glue("{group}2")), as.character) %>%
    distinct_all() %>%
    mutate_at(vars(group, str_glue("{group}2")), factor) %>%
    filter(!!as.name(group) != !!as.name(str_glue("{group}2"))) %>%
    mutate(total = rowSums(select(., .dots = any_of(colS))))

  newLinks <- map(colS %>% c("total"), function(c) {
    x <- quantile(linksMod[[c]], cutoff)
    linksTmp <- linksMod %>%
      select(
        from = all_of(group),
        to = all_of(str_glue("{group}2")),
        weight = !!as.name(c)
      ) %>%
      filter(weight >= x) %>%
      arrange(gtools::mixedsort(from))
  }) %>% setNames(colS %>% c("total"))

  newNodes <- df %>%
    select(any_of(contains(c("camp", "band")))) %>%
    distinct_all() %>%
    rowid_to_column() %>%
    pivot_longer(-rowid) %>%
    mutate_at(vars(name), list(~ str_remove_all(.x, "2"))) %>%
    distinct_all() %>%
    pivot_wider(
      names_from = name,
      values_from = value,
      values_fn = list
    ) %>%
    unchop(c(band, camp)) %>%
    select(-rowid) %>%
    distinct_all() %>%
    select(camp, band) %>%
    mutate_all(factor) %>%
    mutate(name = !!as.name(group))

  return(list(links = newLinks, nodes = newNodes))
}

plotGraph <- function(links,
                      nodes,
                      title = "",
                      weighted = T,
                      color = NA,
                      shape = NA,
                      labels = T,
                      nodeSize = 10,
                      save = F,
                      layout = NA,
                      plotDevice = "svg",
                      groupByBand = F) {
  # layout are coordinates for each node in a data frame format

  # old <- .Random.seed
  # on.exit({
  #   .Random.seed <<- old
  # })
  # set.seed(5010)

  if (is.data.frame(layout)) {
    if (is.matrix(links)) {
      links <- links %>%
        as.data.frame() %>%
        rownames_to_column("from") %>%
        pivot_longer(-from, names_to = "to", values_to = "weight")
    }
    plotdf <- links %>%
      rename(name = from) %>%
      full_join(layout, by = "name") %>%
      left_join(
        layout %>% rename(
          to = name,
          xend = x,
          yend = y
        ),
        by = "to"
      )
  } else if (is.data.frame(links)) {
    n <- igraph::graph_from_data_frame(links, directed = F)
  } else if (is.matrix(links)) {
    colnames(links) <- rownames(links) <- nodes$name
    n <- igraph::graph_from_adjacency_matrix(links, mode = "undirected", weighted = T, diag = F)
  } else {
    stop("check class of links")
  }

  if (!is.data.frame(layout)) {
    if (groupByBand) {
      V(n)$band <- nodes$band
      nG <- contract(n, as.integer(V(n)$band) + 1, vertex.attr.comb = "first")
      V(nG)$name <- nodes$band %>% unique()
      nG <- simplify(nG, edge.attr.comb = "sum")
      plotdf <- ggnetwork::ggnetwork(nG) %>%
        as_tibble()
    } else {
      plotdf <- ggnetwork::ggnetwork(n) %>%
        as_tibble()
    }
  }

  plotdf %<>%
    mutate(weight = weight %>% as.numeric()) %>%
    left_join(nodes %>% mutate_at(vars(name), as.character),
      by = "name"
    )

  gg <-
    plotdf %>%
    ggplot(aes(
      x = x,
      y = y,
      xend = xend,
      yend = yend
    ))
  if (weighted == T) {
    print("adding weights")
    gg <- gg +
      geom_edges(aes(color = weight), size = .75) +
      scale_color_viridis_c(
        direction = -1,
        breaks = c(
          min(plotdf$weight, na.rm = T),
          max(plotdf$weight, na.rm = T)
        ),
        labels = c("low", "high")
      )
  } else {
    gg <- gg +
      geom_edges(
        color = "black",
        size = .5,
        alpha = .5
      )
  }
  if (!is.na(shape) && shape == "band" & is.na(color)) {
    print("adding shape")
    gg <- gg +
      geom_point(aes(shape = !!as.name(shape)),
        size = nodeSize
      ) +
      scale_shape_manual(values = c(15, 17, 19))
  } else if (isTRUE(shape == "band" & !is.na(color))) {
    print("adding shape and color")
    gg <- gg +
      geom_point(
        aes(
          shape = !!as.name(shape),
          fill = !!as.name(color)
        ),
        size = nodeSize
      ) +
      scale_fill_brewer(palette = "Set1") +
      scale_shape_manual(values = c(21, 22, 24))
  } else if (shape != "band" & !is.na(shape)) {
    stop("too many shapes -- maybe try color?")
  } else if (is.na(shape) & !is.na(color)) {
    print("adding color")
    gg <- gg +
      geom_point(aes(fill = !!as.name(color)),
        size = nodeSize,
        shape = 21
      ) +
      scale_fill_brewer(palette = "Set1")
  } else {
    gg <- gg +
      geom_point(
        size = nodeSize,
        color = "black",
        shape = 19
      )
  }

  if (labels == T) {
    print("adding labels")
    gg <- gg +
      geom_nodetext(aes(label = name),
        # size = 3,
        fontface = "bold",
        color = "white"
      )
  }
  # determine if legend shape should be overridden
  if (!is.na(shape)) {
    f1 <- guide_legend()
  } else if (isTRUE(!is.na(color) && shape != color) | !is.na(color)) {
    f1 <- guide_legend(override.aes = list(shape = 21))
  } else {
    f1 <- guide_legend()
  }
  gg <- gg +
    theme_blank() +
    theme(
      legend.position = "bottom",
      legend.background = element_rect(fill = F),
      plot.margin = margin(0, 0, 0, 0)
    ) +
    coord_cartesian(clip = "off") +
    guides(
      fill = f1,
      shape = guide_legend(direction = "horizontal"),
      color = guide_colorbar(
        direction = "horizontal",
        title.position = "top"
      )
    ) +
    ggtitle(title)
  print(gg)
  if (save == T) {
    if (plotDevice == tolower("svg")) {
      ggsave(
        filename = glue::glue("plots/{title}.svg"),
        plot = gg
      )
    } else if (plotDevice == tolower("png")) {
      ggsave(
        filename = glue::glue("plots/{title}.png"),
        plot = gg,
        dpi = 450
      )
    } else {
      stop("how about svg or png instead?")
    }
  }
  return(gg)
}

#' Calculate weighted Jaccard
weightedJaccard <- function(m) {
  #' weighted jaccard similarity matrix setup
  #' from:
  #' https://rpubs.com/lgadar/weighted-jaccard

  r <- matrix(0, nrow = nrow(m), ncol = nrow(m))
  rownames(r) <- rownames(m)
  colnames(r) <- rownames(m)

  # weighted jaccard
  pairs <- t(combn(1:nrow(m), 2))
  for (i in 1:nrow(pairs)) {
    num <-
      sum(sapply(1:ncol(m), function(x) {
        (min(m[pairs[i, 1], x], m[pairs[i, 2], x]))
      }))
    den <-
      sum(sapply(1:ncol(m), function(x) {
        (max(m[pairs[i, 1], x], m[pairs[i, 2], x]))
      }))
    r[pairs[i, 1], pairs[i, 2]] <- num / den
    r[pairs[i, 2], pairs[i, 1]] <- num / den
  }
  r[which(is.na(r))] <- 0
  diag(r) <- 0
  return(r)
}

#' Format netlogo adjacency tables as matrices
#'
#' @param x Netlogo matrix
#'
#' @return
#' @export
#'
#' @examples
string2MatrixR <- function(x) {
  m <- x %>%
    str_remove_all("\\[\\[|\\]\\]") %>%
    str_split_1(pattern = fixed("] [")) %>%
    as.matrix() %>%
    apply(1, str_split_1, pattern = " ")
  m[] <- as.double(m[])
  return(m)
}

library(Rcpp)
sourceCpp("analysis/string2matrix.cpp")

#' Format Netlogo Adjacency Tables
#'
#' @param p run number
#' @param nB number of bands
#' @param cutoff cutoff value
#' @param x adjacency table formatted as string
#'
#' @return
#' @export
#'
#' @examples
formatNLAdjTables <- function(x, p = NULL, standardize = T) {
  result <- purrr::map(x, function(x1) {
    m <- string2matrix(x1)
    nB <- ncol(m) / 3
    nms <- nB:(nB * 3 + nB - 1)
    bands <- data.frame(band = as.character(sort(rep(0:(nB - 1), 3))), camp = as.character(nms))
    colnames(m) <- rownames(m) <- nms
    df <- as.data.frame(m) %>%
      rownames_to_column("from") %>%
      pivot_longer(all_of(nms %>% as.character()),
        names_to = "to",
        values_to = "weight"
      )
    if (!is.null(p)) {
      df <- df %>%
        mutate(property = p) %>%
        mutate(id = map2_chr(from, to, function(f, t) {
          sort(c(f, t)) %>%
            paste(collapse = "_")
        })) %>%
        distinct(id, .keep_all = T) %>%
        select(-id) %>%
        separate(
          col = property,
          into = c("run_number", "property"),
          sep = "_",
          extra = "drop"
        )
    }

    df <- df %>%
      rowwise() %>%
      mutate(id = sort(c(from, to)) %>% paste(collapse = "_")) %>%
      ungroup() %>%
      distinct(id, .keep_all = T) %>%
      filter(from != to) %>%
      filter(weight > 0)

    if (standardize) {
      df <- df %>%
        mutate(weight = weight / max(weight))
    }
    df <- df %>%
      # mutate_at(vars(from,to),as.character) %>%
      left_join(bands %>% rename(from = camp), by = "from")
    return(df)
  })

  # nodes = tibble(camp = nms,
  #                band = rep((1:nB - 1), 3) %>% sort) %>%
  #   mutate_all(as.character) %>%
  #   mutate(name = camp)
  # return(list(links = links, nodes = nodes))
  return(result)
}

#' Format assemblages into a data frame structure
formatAssemblages <- function(assemblages) {
  # tictoc::tic()
  tmp <- assemblages %>%
    # remove the aggregation sites
    filter(!is.na(band)) %>%
    mutate_at(vars(band, site), as.character)
  pot <- tmp %>%
    select(band, site, potAssemblage) %>%
    mutate(potAssemblage = potAssemblage %>%
      map(bracket2list)) %>%
    unnest(potAssemblage) %>%
    mutate(type = "pot") %>%
    select(band, site, tick, type, tool)
  point <- tmp %>%
    select(band, site, pointAssemblage) %>%
    mutate(pointAssemblage = pointAssemblage %>%
      map(bracket2list)) %>%
    unnest(pointAssemblage) %>%
    mutate(type = "point") %>%
    select(band, site, tick, type, tool)
  # print(tictoc::toc()$elapsed)
  return(bind_rows(pot, point))
}

#' Create attribute network
attrNetworks <- function(runResults,
                         attr = "pointAttributes",
                         group = "who") {
  links <- runResults$nodes %>%
    select(who, band, camp, !!as.name(attr)) %>%
    mutate(attr = !!as.name(attr) %>% str_remove_all("\\[|\\]")) %>%
    select(-!!as.name(attr)) %>%
    separate_rows(attr, sep = " ")
  # table(links$who,links$attr)
  if (group == "camp") {
    links %<>%
      select(-who) %>%
      distinct_all()
  } else if (group == "band") {
    links %<>%
      select(-who, -camp) %>%
      distinct_all()
  }
  d <-
    table(
      links[[group]],
      links$attr
    )
  n <- nrow(d)
  m <- matrix(nrow = n, ncol = n)
  rownames(m) <- rownames(d)
  colnames(m) <- rownames(d)
  for (i in 1:n) {
    for (j in 1:n) {
      r1 <- which(d[i, ] > 0)
      r2 <- which(d[j, ] > 0)
      m[i, j] <- intersect(
        r1,
        r2
      ) %>%
        length()
    }
  }
  m
  diag(m) <- 0
  m[which(m < 2)] <- 0
  g <- list()
  g$links <- as.data.frame(m) %>%
    rownames_to_column("from") %>%
    pivot_longer(2:ncol(.), names_to = "to", values_to = "weight") %>%
    filter(weight > 0) %>%
    mutate(id = map2_chr(
      from, to,
      ~ paste0(min(.x, .y), "_", max(.x, .y))
    )) %>%
    distinct(id, .keep_all = T) %>%
    select(-id)
  g$nodes <- links %>%
    mutate(name = !!as.name(group)) %>%
    mutate_all(as.character) %>%
    mutate_at(vars(any_of(c(
      "who", "name", "camp", "band"
    ))), list(function(x) {
      factor(x, levels = unique(x) %>% gtools::mixedsort())
    }))
  return(g)
}

#' Select a proportion of the data either randomly or nonrandomly
sampleData <- function(df,
                       p, # proportion of each site to keep
                       method # c("random","nonrandom")
) {
  s <- df %>%
    group_by(site) %>%
    count() %>%
    mutate(
      c = floor(n * p),
      c = case_when(c == 0 ~ 1, TRUE ~ c)
    )
  if (method == "random") {
    rs <- df %>%
      group_by(site) %>%
      group_split()
    r <- map2_dfr(rs, s$c, ~ {
      .x %>%
        slice_sample(n = .y)
    })
  } else {
    rs <- df %>%
      slice(gtools::mixedorder(site)) %>%
      arrange(site) %>%
      rowid_to_column() %>%
      group_by(site) %>%
      group_split()
    r <- map2_dfr(
      rs,
      s$c, function(rsi, si) {
        print("asdf")
        indx <- rsi %>%
          group_by(site) %>%
          slice_sample(n = 1) %>%
          pull(rowid)

        indx2 <- indx + si
        if (indx2 > max(rsi$rowid)) {
          tmp1 <- rsi %>%
            filter(
              rowid >= indx,
              rowid <= max(rsi$rowid)
            )
          n <- si - nrow(tmp1)
          indx <- rsi$rowid[which.min(rsi$rowid)]
          indx2 <- indx + n

          tmp2 <- rsi %>%
            filter(
              rowid >= indx,
              rowid < indx2
            )
          rsi <- bind_rows(tmp1, tmp2)
        } else {
          rsi %<>%
            filter(
              rowid >= indx,
              rowid < indx2
            )
        }
        return(rsi)
      }
    )
  }
  return(r)
}

compareMissing <- function(df,
                           group,
                           prop, # proportion of each site to keep
                           method # c("random","nonrandom")){
) {
  df %<>% rename(group = !!as.name(group))
  original <- table(df$site, df$tool) %>%
    prop.table(1) %>%
    {
      . * 100
    } %>%
    weightedJaccard()
  results <- map_dfr(prop, function(p) {
    print(paste("running", method, p))
    reduced <- sampleData(df = df, p = p, method = method)
    sampled <- table(reduced$site, reduced$tool) %>%
      prop.table(1) %>%
      {
        . * 100
      } %>%
      weightedJaccard()
    l <- list(original, sampled)
    distance <- try(round(nd.dsd(l,
      out.dist = T,
      type = "Lap"
    )$D, 3))
    if (inherits(distance, "try-error")) {
      distance <- NA_real_
    }
    result <- tibble(
      count = nrow(reduced),
      proportion = p,
      method = method,
      distance = unlist(distance)[1]
    )
    return(result)
  })
  return(results)
}

formatMatrix <- function(net, q = 0, binarize = F, rescale = T, groupByBand = F) {
  if (inherits(net, "list")) {
    net <- net[[1]]
  }
  if (inherits(net[1, 1], "character")) {
    net <- apply(net, 1, as.numeric)
  }
  if (groupByBand) {
    nB <- ncol(net) / 3
    bNms <- 0:(nB - 1)
    nms <- (nB):(ncol(net) + (nB - 1))
    colnames(net) <- rownames(net) <- nms
    n <- igraph::graph_from_adjacency_matrix(net, mode = "undirected", weighted = T, diag = F)
    igraph::V(n)$band <- bNms
    nG <- igraph::contract(n, as.integer(V(n)$band) + 1, vertex.attr.comb = "first")
    igraph::V(nG)$name <- bNms
    nG <- igraph::simplify(nG, edge.attr.comb = "sum")
    net <- igraph::as_adjacency_matrix(nG, attr = "weight") %>% as.matrix()
  }
  if (q > 0) {
    cutoff <- quantile(net, q)
    net[which(net <= cutoff)] <- 0
    if (binarize) {
      net[which(net > cutoff)] <- 1
    }
  }
  if (rescale) {
    net <- log(net + 1)
    net <- scales::rescale(net)
  }
  rownames(net) <- colnames(net) <- 1:ncol(net)
  return(net)
}

adjacency_compare <- function(net1, net2, q = 0, binarize = F, rescale = T) {
  net1 <- formatMatrix(net1, q = q, binarize = binarize, rescale = rescale)
  net2 <- formatMatrix(net2, q = q, binarize = binarize, rescale = rescale)
  net_diff <- abs(net1 - net2)
  out <-
    1 - (sum(net_diff) / (nrow(net_diff) * (nrow(net_diff) - 1)))
  return(out)
}

adjacency_compare2 <- function(net1, net2) {
  if (!identical(dim(net1), dim(net2))) {
    warning("matrices must be the same size: returning NA")
    return(NA)
  } else if (any(is.na(net1)) || any(is.na(net2))) {
    warning("matrices must not be NA: returning NA")
    return(NA)
  } else if (any(is.infinite(net1)) || any(is.infinite(net2))) {
    warning("matrices must not be Infinite: returning NA")
    return(NA)
  } else {
    diag(net1) <- diag(net2) <- 0
    net_diff <- abs(net1 - net2)
    out <- 1 - sum(net_diff) / (nrow(net_diff) * nrow(net_diff))
    return(out)
  }
}

adjacency_compare3 <- function(net1, net2) {
  net_diff <- sum(abs(net1 - net2))
  return(net_diff)
}

cor_compare <- function(net1, net2) {
  tryCatch(
    {
      if (!identical(dim(net1), dim(net2))) {
        warning("matrices must be the same size: returning NA")
        return(NA)
      } else if (any(is.na(net1)) || any(is.na(net2))) {
        warning("matrices must not be NA: returning NA")
        return(NA)
      } else if (any(is.infinite(net1)) || any(is.infinite(net2))) {
        warning("matrices must not be Infinite: returning NA")
        return(NA)
      } else {
        diag(net1) <- diag(net2) <- 0
        out <- sna::gcor(net1, net2)
        return(out)
      }
    },
    error = function(e) {
      return(NA)
    }
  )
}

qapTestFun <- function(net1, net2) {
  tryCatch(
    {
      if (!identical(dim(net1), dim(net2))) {
        warning("matrices must be the same size: returning NA")
        return(NA)
      } else if (any(is.na(net1)) || any(is.na(net2))) {
        warning("matrices must not be NA: returning NA")
        return(NA)
      } else if (any(is.infinite(net1)) || any(is.infinite(net2))) {
        warning("matrices must not be Infinite: returning NA")
        return(NA)
      } else {
        diag(net1) <- diag(net2) <- 0
        out <- sna::qaptest(list(net1, net2), sna::gcor, g1 = 1, g2 = 2)$pleeq
        return(out)
      }
    },
    error = function(e) {
      return(NA)
    }
  )
}

adjacency_correlation <- function(net1, net2, q = 0, binarize = T, rescale = F) {
  net1 <- formatMatrix(net1, q = q, binarize = binarize, rescale = rescale)
  net2 <- formatMatrix(net2, q = q, binarize = binarize, rescale = rescale)
  library(igraph)
  g1 <- graph_from_adjacency_matrix(net1, mode = "undirected", diag = F)
  g2 <- graph_from_adjacency_matrix(net2, mode = "undirected", diag = F)
  library(multinet)
  net <- ml_empty("comparison")
  add_igraph_layer_ml(net, g1, "net1")
  add_igraph_layer_ml(net, g2, "net2")
  m <- layer_comparison_ml(net, method = "pearson.degree", mode = "all")
  diag(m) <- NA
  return(max(m, na.rm = T))
}

adjacency_laplacian <- function(net1, net2, q = 0, binarize = F, rescale = T) {
  net1 <- formatMatrix(net1, q = q, binarize = binarize, rescale = rescale)
  net2 <- formatMatrix(net2, q = q, binarize = binarize, rescale = rescale)
  library(NetworkDistance)
  l <- list(net1, net2)
  m <- round(nd.dsd(l,
    out.dist = T,
    type = "Lap"
  )$D, 3)
  return(m)
}

adjacency_laplacian2 <- function(net1, net2) {
  if (!identical(dim(net1), dim(net2))) {
    warning("matrices must be the same size: returning NA")
    return(NA)
  } else if (any(is.na(net1)) || any(is.na(net2))) {
    warning("matrices must not be NA: returning NA")
    return(NA)
  } else if (any(is.infinite(net1)) || any(is.infinite(net2))) {
    warning("matrices must not be Infinite: returning NA")
    return(NA)
  } else {
    library(NetworkDistance)
    l <- list(net1, net2)
    m <- round(nd.dsd(l,
      out.dist = T,
      type = "Lap"
    )$D, 3)
    return(m)
  }
}


compareAdjacencyResults <- function(runs) {
  results <- runs %>%
    mutate_all(as.character) %>%
    pivot_longer(-run_number,
      names_to = "property",
      values_to = "value"
    ) %>%
    mutate(property = str_remove_all(
      property,
      "matrix_to_row_list_camp_|_adjacency"
    )) %>%
    pivot_wider(names_from = property, values_from = value) %>%
    pivot_longer(visiting:pot, names_to = "property", values_to = "adj") %>%
    unite("id", run_number, property, remove = F)

  nB <- runs %>%
    select(n_of_bands) %>%
    dplyr::pull() %>%
    unique()

  formatted <- map2_df(
    results$adj, results$id,
    formatNLAdjTables, nB
  )
  nms <- nB:(nB * 3 + nB - 1)
  nodes <- tibble(
    camp = nms,
    band = rep((1:nB - 1), 3) %>% sort()
  ) %>%
    mutate_all(as.character) %>%
    mutate(name = camp)

  groups <- formatted$property %>%
    unique() %>%
    sort()
  comparison <- tibble()
  for (i in formatted$run_number %>%
    unique() %>%
    as.integer()) {
    print(i)
    try({
      run <- formatted %>%
        filter(run_number == i)
      nms <- run$property %>%
        unique() %>%
        sort()
      run %<>%
        group_by(property) %>%
        group_split() %>%
        setNames(nms)


      # plotGraph(run$hunting,nodes, "hunting",shape = 'band', color = 'band')
      cutoff <- quantile(run$learning$weight, 0.5)
      run$learning %<>% filter(weight > cutoff)
      # plotGraph(run$learning,nodes, "learning",shape = 'band', color = 'band')
      cutoff <- quantile(run$trading$weight, 0.5)
      run$trading %<>% filter(weight > cutoff)
      # plotGraph(run$trading,nodes, "trading",shape = 'band', color = 'band')
      cutoff <- quantile(run$visiting$weight, 0.5)
      run$visiting %<>% filter(weight > cutoff)
      # plotGraph(run$visiting,nodes, "visiting",shape = 'band', color = 'band')
      cutoff <- quantile(total$weight, 0.5)
      # plotGraph(total,nodes, "total",shape = 'band', color = 'band')


      cutoff <- quantile(run$point$weight, 0.5)
      run$point %<>% filter(weight > cutoff)
      # plotGraph(run$point,nodes, "point",shape = 'band', color = 'band')
      cutoff <- quantile(run$pot$weight, 0.5)
      run$pot %<>% filter(weight > cutoff)
      # plotGraph(run$pot,nodes, "pot",shape = 'band', color = 'band')

      # multilayer network
      library(igraph)
      library(multinet)
      net <- ml_empty("comparison")
      ig1 <- graph_from_data_frame(d = run$pot, directed = F, vertices = nodes)
      E(ig1)$weight <- run$pot$weight
      add_igraph_layer_ml(n = net, g = ig1, name = "pots")
      ig2 <- graph_from_data_frame(d = run$point, directed = F, vertices = nodes)
      E(ig2)$weight <- run$point$weight
      add_igraph_layer_ml(n = net, g = ig2, name = "points")
      ig3 <- graph_from_data_frame(d = run$hunting, directed = F, vertices = nodes)
      E(ig3)$weight <- run$hunting$weight
      add_igraph_layer_ml(n = net, g = ig3, name = "hunting")
      ig4 <- graph_from_data_frame(d = run$learning, directed = F, vertices = nodes)
      E(ig4)$weight <- run$learning$weight
      add_igraph_layer_ml(n = net, g = ig4, name = "learning")
      ig5 <- graph_from_data_frame(d = run$trading, directed = F, vertices = nodes)
      E(ig5)$weight <- run$trading$weight
      add_igraph_layer_ml(n = net, g = ig5, name = "trading")
      ig6 <- graph_from_data_frame(d = run$visiting, directed = F, vertices = nodes)
      E(ig6)$weight <- run$visiting$weight
      add_igraph_layer_ml(n = net, g = ig6, name = "visiting")
      m2raw <- layer_comparison_ml(net, method = "pearson.degree", mode = "all")
      m2 <- round(m2raw, 2) %>%
        as.data.frame() %>%
        rownames_to_column("from") %>%
        pivot_longer(-from, names_to = "to") %>%
        filter(from != to) %>%
        rowwise() %>%
        mutate(variable = paste0(
          paste(sort(c(from, to)), collapse = "->"),
          "-pearson.degree"
        ), .before = 1) %>%
        select(-from, -to) %>%
        distinct_all()
      # m %>% knitr::kable()

      ig1m <- igraph::as_adjacency_matrix(ig1)
      ig2m <- igraph::as_adjacency_matrix(ig2)
      ig3m <- igraph::as_adjacency_matrix(ig3)
      ig4m <- igraph::as_adjacency_matrix(ig4)
      ig5m <- igraph::as_adjacency_matrix(ig5)
      ig6m <- igraph::as_adjacency_matrix(ig6)

      result <- bind_rows(
        tibble(
          variable = "pots->points-AdjComparison",
          value = adjacency_compare(ig1m, ig2m)
        ),
        tibble(
          variable = "pots->hunting-AdjComparison",
          value = adjacency_compare(ig1m, ig3m)
        ),
        tibble(
          variable = "pots->learning-AdjComparison",
          value = adjacency_compare(ig1m, ig4m)
        ),
        tibble(
          variable = "pots->trading-AdjComparison",
          value = adjacency_compare(ig1m, ig5m)
        ),
        tibble(
          variable = "pots->visiting-AdjComparison",
          value = adjacency_compare(ig1m, ig6m)
        ),
        tibble(
          variable = "points->hunting-AdjComparison",
          value = adjacency_compare(ig2m, ig3m)
        ),
        tibble(
          variable = "points->learning-AdjComparison",
          value = adjacency_compare(ig2m, ig4m)
        ),
        tibble(
          variable = "points->trading-AdjComparison",
          value = adjacency_compare(ig2m, ig5m)
        ),
        tibble(
          variable = "points->visiting-AdjComparison",
          value = adjacency_compare(ig2m, ig6m)
        ),
        tibble(
          variable = "hunting->learning-AdjComparison",
          value = adjacency_compare(ig3m, ig4m)
        ),
        tibble(
          variable = "hunting->trading-AdjComparison",
          value = adjacency_compare(ig3m, ig5m)
        ),
        tibble(
          variable = "hunting->visiting-AdjComparison",
          value = adjacency_compare(ig3m, ig6m)
        ),
        tibble(
          variable = "learning->trading-AdjComparison",
          value = adjacency_compare(ig4m, ig5m)
        ),
        tibble(
          variable = "learning->visiting-AdjComparison",
          value = adjacency_compare(ig4m, ig6m)
        ),
        tibble(
          variable = "trading->visiting-AdjComparison",
          value = adjacency_compare(ig5m, ig6m)
        ),
        m2
      )
      vars <- results %>%
        filter(run_number == i) %>%
        select(-id, -property, -adj) %>%
        distinct_all() %>%
        pivot_longer(-run_number, names_to = "variable")
      result <- vars %>%
        bind_rows(result %>%
          mutate_at(vars(value), round, digits = 2) %>%
          mutate(run_number = i) %>%
          mutate_all(as.character))
      comparison <- bind_rows(comparison, result)
    })
  }
  return(comparison)
}

matrix2rank <- function(adj) {
  m <- string2matrix(adj)
  v <- sort(unique(as.vector(m)), decreasing = T)
  vRank <- 1:length(v)
  m[] <- match(m, v)
  diag(m) <- 0
  return(m)
}

# net1 = matrix2rank(runs$wordMatrixToRowListCampLearningAdjacency[2])
# net2 = matrix2rank(runs$wordMatrixToRowListCampPotAdjacency[2])
# adjacency_compare2(net1,net2)

matrix2rank2 <- function(adj) {
  m <- adj
  v <- sort(unique(as.vector(m)), decreasing = T)
  vRank <- 1:length(v)
  m[] <- match(m, v)
  diag(m) <- 0
  return(m)
}

matrix2zscore <- function(adj) {
  m <- string2matrix(adj)
  v <- (m - mean(m)) / sd(m)
  m[] <- v
  diag(m) <- 0
  m[] <- round(m, 3)
  m <- scales::rescale(m)
  return(m)
}

matrix2zscoreBand <- function(adj) {
  m <- string2matrix(x)
  nc <- ncol(m) / 3
  bands <- rep(0:(nc - 1), 3) %>% sort()
  camps <- 1:length(bands)
  dict <- tibble(bands, camps)
  m <- as.data.frame(m) %>%
    rowid_to_column("from") %>%
    pivot_longer(-from, names_to = "to") %>%
    mutate(to = str_remove_all(to, "V") %>% as.integer()) %>%
    left_join(dict, by = c("from" = "camps")) %>%
    rename(bands1 = bands) %>%
    left_join(dict, by = c("to" = "camps")) %>%
    rename(bands2 = bands) %>%
    group_by(bands1, bands2) %>%
    summarize(value = sum(value)) %>%
    ungroup() %>%
    pivot_wider(names_from = bands2, values_from = value) %>%
    column_to_rownames("bands1") %>%
    as.matrix()
  diag(m) <- 0
  v <- (m - mean(m)) / sd(m)
  m[] <- v
  diag(m) <- 0
  m[] <- round(m, 3)
  m <- scales::rescale(m)
  return(m)
}

matrix2one <- function(adj) {
  m <- string2matrix(adj)
  m[] <- m / max(m)
  diag(m) <- 0
  m[] <- round(m, 3)
  return(m)
}

# net1 = matrix2rank(runs$wordMatrixToRowListCampLearningAdjacency[2])
# net2 = matrix2rank(runs$wordMatrixToRowListCampPotAdjacency[2])
# adjacency_compare2(net1,net2)

formatNLAdjTablesOld <- function(x, p, nB = 3) {
  nms <- nB:(nB * 3 + nB - 1)
  links <- x %>%
    strsplit(" |\\[|\\]") %>%
    unlist() %>%
    as.numeric() %>%
    .[which(!is.na(.))] %>%
    matrix(
      nrow = length(nms),
      ncol = length(nms),
      byrow = T
    ) %>%
    set_colnames(nms) %>%
    set_rownames(nms) %>%
    as.data.frame() %>%
    rownames_to_column("from") %>%
    pivot_longer(all_of(nms %>% as.character()),
      names_to = "to",
      values_to = "weight"
    ) %>%
    mutate(
      property = p,
      id = map2_chr(from, to, function(f, t) {
        sort(c(f, t)) %>%
          paste(collapse = "_")
      })
    ) %>%
    distinct(id, .keep_all = T) %>%
    select(-id) %>%
    separate(
      col = property,
      into = c("run_number", "property"),
      sep = "_",
      extra = "drop"
    ) %>%
    filter(from != to) %>%
    filter(weight > 0)
  return(links)
}
#' Format assemblages into a data frame structure

eigenDiff <- function(adj1, adj2) {
  m1 <- matrix2rank(adj1)
  m2 <- matrix2rank(adj2)
  e1 <- sna::evcent(m1)
  e2 <- sna::evcent(m2)
  return(sum(abs(e1 - e2)))
}

# tmp = runs %>% slice(1:100)

# subset common vertices from 2 networks
network_subset_common <-
  function(net1, net2, namevar = "vertex.names") {
    net1_names <-
      network::get.vertex.attribute(net1, attrname = namevar)
    net2_names <-
      network::get.vertex.attribute(net2, attrname = namevar)
    common_names <- net1_names[which(net1_names %in% net2_names)]
    out1 <- net1 %s% which(net1 %v% namevar %in% common_names)
    out2 <- net2 %s% which(net2 %v% namevar %in% common_names)
    return(list(out1, out2))
  }

R Environment

R version=4.4.3 callr=3.7.6 ergm=4.8.1 network=1.19.0 mgcv=1.9-1 nlme=3.1-167 Rcpp=1.0.14 ggthemes=5.1.0 magrittr=2.0.3 BayesFactor=0.9.12-4.7 Matrix=1.7-2 coda=0.19-4.1 magick=2.8.5 tictoc=1.2.1 ggnetwork=0.5.13 igraph=2.1.4 ggpubr=0.6.0 arrow=19.0.1 lubridate=1.9.4 forcats=1.0.0 stringr=1.5.1 purrr=1.0.4 readr=2.1.5 tidyr=1.3.1 tibble=3.2.1 ggplot2=3.5.1 tidyverse=2.0.0 here=1.0.1 flextable=0.9.7 gt=0.11.1 dplyr=1.1.4 rio=1.2.3

ODD

“ArchMatNet: Archaeological Material Networks”

ODD: Overview, Design concepts, and Details

Cecilia Padilla-Iglesias, Claudine Gravel-Miguel, and Robert Bischoff1

23 December 2022

Purpose

The purpose of the model is to investigate how different factors affect the ability of researchers to reconstruct prehistoric social networks from artifact stylistic similarities, as well as the overall diversity of cultural traits observed in archaeological assemblages. Given that cultural transmission and evolution are affected by multiple interacting phenomena, our model allows to simultaneously explore six sets of factors that may condition how social networks relate to shared culture between individuals and groups:

  1. Factors relating to the structure of social groups
  2. Factors relating to the cultural traits in question
  3. Factors relating to individual learning strategies
  4. Factors relating to the environment
  5. Factors relating to the context in which different types of cultural traits are learned
  6. Factors relating to the method used to reconstruct ancient social networks

Process overview and scheduling

At the beginning of each time-step, the model checks what type of environment is set. If it is only one large biome, it rolls the die to see if the environment will change. If it does, the new environment is a random value between 0 and 99. Patches get that new environment.

If the user wants highResolutionData, the model checks if it is time to export data (every 1000 time-steps). If it is the case, it calls procedures that compile data and the procedure that exports it to the correct CSV file.

People remove any superfluous alliances they may have accumulated in previous tick, using the cleanUpAlliances procedure.

People run the possibility to invent new types and new variants and create one pot and point each. A person may create a new pot and/or point each tick. Each new object is created based on a randomly selected trait from the person’s trait list. Each object also records the tick when it was created for use in subsequent analyses.

The model then checks if it is time to aggregate (based on the aggregation-freq slider). If it is, it sets aggregation to T, and calls the aggregate procedure, which moves people to their respective aggregation camps.

While aggregating, the model keeps track of the aggregation’s length. If it has reached the aggregationDuration limit, people go back to their own camp. When they are aggregating, people “visit” with all the people they are aggregated with, which is recorded in the interactions’ “visiting” counter, as well as the inter-camps campInteractions “visiting” counters. Aggregated people can go hunting with a certain probability.

When an aggregation is not occurring, people can do different activities, with a given probability. They can create a new alliance with someone from another camp. They can go visit an ally. They can go hunting.

When hunting, people gather other people from their band and go to a different location together. When the hunt occurs during an aggregation, the hunting site is within 15 patches of the aggregation site. When the hunt occurs at other times, the hunting site is simply within the band territory. On a hunt, the hunters calculate their pointPrestige, and have different probabilities to drop a point, trade points with one another, and transmit trait knowledge to one another. Hunters keep track of how many times they hunted together. This interaction is also recorded by the campInteractions counter between camps. When the hunt is over, everyone returns to their own camp.

When visiting, one person chooses one ally to go visit at random. They move to that ally’s camp and visit with all the allies present there. The people who are involved in the visit keep track of how many times they visited one another, which is also recorded by the campInteractions visiting counter. The people who participate in the visit have different probabilities to drop one pot, drop one point, trade objects with one another, and transmit trait knowledge to one another. When those activities are done, the visitor returns to their own camp.

If recordPrestige is T, camps call the recordCampPrestige procedure. This procedure helps to follow the progression of mean and max prestige of the people of each camp, to see how those values change over time.

All people have a certain probability to call the loseTrait procedure, in which they forget one pot or point trait.

All people call the transmitCulture procedure to share knowledge of traits, and tradeItems to trade objects from their inventories.

Then, all people call the dropItem procedure, in which they have a certain probability to drop one pot or one point.

After the first 1000 ticks, people check if they want to migrate or not at every tick, based on a probability set by the user. When genderedMigration is T, only people with gender 1 roll the die and they have twice the probability to migrate as everyone does when genderedMigration is F. When migration occurs, the person who is migrating chooses another person to switch place with. When genderedMigration is T, only people of gender 1 can be selected. The two people switch home camps.

The model checks to see if it has reached the maximum number of ticks. If it has and the user wants to see the metrics, it computes and prints all the metrics to the Command Center, exports results as a CSV file if saveResults is T, and stops.

Design concepts

Basic principles

We built this model to allow testing assumptions that archaeologists have made about the possibility to reconstruct social networks from material culture, as well as improve our methods to reconstruct said networks.

Emergence

We expect styles to emerge over time within well-connected people who share the same traits and innovate on the shared known traits.

Adaptation

The people do not adapt to their environment or their context. Their behavior is constant throughout the run.

Objectives

One objective people have is to learn fitter traits (when learning via Prestige).

Learning

People regularly learn from the people around them. This is the base of the cultural transmission subprocedure. While the learning is not considered adaptation, people will learn the traits that are most popular at any given time, traits which change as innovations gain popularity or are lost.

Prediction

There is no prediction.

Sensing

People sense the environment onto which they hunt, and it affects their prestige. When meeting one another, people can sense which traits the others know. People also know who their allies are, and how often they have visited them. Camps and bands sense the knowledge of their inhabitants. None of the knowledge involve errors. Interaction Interactions make the most important part of the model. People hunt together and visit one another. When interacting, people can trade items or share knowledge of some traits.

Stochasticity

To keep the model general, we use stochasticity often. A lot of the activities are done with a certain probability, controlled by Interface sliders. Traded items are chosen at random amongst traders’ inventories. When creating a new variant, the type to innovate upon is chosen randomly amongst the traits known by the innovator. Forgotten traits are also selected at random. This use of stochasticity is also because we do not know how often and how those activities were done in the past. Therefore, it is simpler to use probabilities to get a wider range of behaviors.

Collectives

People are grouped by camps, which are grouped by bands. There are three camps per band. People aggregate and hunt with the people from their band, therefore, their group affect who they will interact with the most, which in turn affects the similarity of their known traits and created objects.

Observation

The model can produce different types of outputs. It can calculate and print the following metrics (in the following names, X can be point or pot): 1) AssemblageSizeX: the sum of pots or points discarded at camps and bands; 2) interactionFreqY (where Y can be hunting, visiting, learning, and trading): the sum of all Y interactions between people; 3) brDistBandX: mean Brainerd Robinson distance between known traits, grouped by bands; 4) brDistCampX: mean Brainerd Robinson distance between known traits, grouped by camps; 5) meanPeoplePrestige: mean hunting prestige of all people; 6) modBandX: modularity metric which evaluates if there is more between-band links than within-band links, based on the Jaccard similarity of traits known by people; 7) modCampX: modularity metric which evaluates if there is more between-camp links than within-camp links, based on the Jaccard similarity of traits known by people; 8) networkMetricX: difference between the clustering coefficient of the network created by people who interacted with one another (hunting and learning) and the clustering coefficient of the network created by similarities in objects (point and pot); 9) PhiSTXBand: PhiST measure to see if variation in known traits is structured around bands; 10) PhiSTXCamp: PhiST measure to see if variation in known traits is structured around camps; 11) nTraitX: sum of the known traits at the end of the simulation (includes duplicates); and 12) elapsedTime: time elapsed since the start of the run.

In addition, the model can produce CSV outputs recording the parameter values, all the interactions between people, as well as the characteristics of all objects discarded at camp/band sites.

Input data

This model does not import external data.

Submodels

Social activities

People meet one another through hunting and visiting. Hunting is done with bandmates only, but visiting is done between allies, which can be of other bands, as hunter-gatherers in the ethnographic literature tend to travel further distances to visit acquaintances than to carry out subsistence activities (Hewlett et al. 1982; Weissner, 1977; Hill et al. 2014)

When someone initiates a hunt, they select a certain number of bandmates to bring with them. They all move to the hunting spot, which is either within the band territory, or within 15 patches of the aggregation site (when the hunt happens during an aggregation). All people taking part in one hunt update their and their camps’ interaction “hunting” counter. Hunting people calculate the fitness of their known traits in the current environment, by calculating how many of their point traits were invented in the current environment. They then have a certain probability of dropping a point, trading items and/or transmit trait knowledge with a fellow hunter. After that, all hunters return to their own camp.

To visit an ally, a person has first to have allies. Allies are formed randomly and with a given probability at each tick. When forming an alliance, a person chooses another person at random. That ally can be of any camp other than ego’s. When someone has more than the limit of allies allowed, they remove one. To select which alliance will go, the person uses weighted sampling, with the number of visits serving as the inverted weight (alliances of people who visit one another often will have less chances of being removed).

When visiting an ally, a person moves to the camp of one of their allies, choosing at random. If they have other allies in that camp, they all participate in the visit and they all update their and their camps’ interaction “visiting” counter. They then have a certain probability of dropping a point and a pot, trading items and/or transmit trait knowledge with a fellow ally. After that, the visitor returns to their own camp. Item exchange

The tradeItems procedure controls exchange between people. It is only run if the person has a point or pot to give and has found a trading partner on the same patch who also has an object to trade in its inventory. This is done first for points, and then for pots. The trading partners for both objects do not need to be the same, but they can be. Both partners choose one item at random from their inventory, and swap them. The interaction link between the agents and the campInteraction links between their camps record that interaction.

Cultural transmission

Cultural transmission happens when the transmitCulture is called (in GO, hunt, and visit). This procedure identifies if the transmission will take place (based on probabilities) and the type of transmission. For pots, and then points, people roll the die to see if they will transmit information. When successful, they will call transmitPrestige with a probability set by the learningMethod slider. Every other time, they will call transmitConformism. These procedures return a set of new traits, following the rules detailed below.

After learning the new traits, people craft a new item (calling the craftItem procedure) and remove any duplicate traits they might have just learned. The transmitConformism [ toolType ] reporter takes as toolType either “pot” or “point” and is used to determine which attributes are copied.

  1. Ego asks all people on the same patch (including itself) to select their traits of the given toolType.
  2. If visibility changes how often traits are transmitted, everyone uses the getMostlyVisibleTraits reporter to get mostly visible traits from their list. High visibility traits are chosen for transmission 98% of the time. This value is arbitrary but generally fits expectations from a broad reading of the ethnographic literature.
  3. Everyone pools their list together and the reporter reports the most common trait. Then, the interaction links between all people update their learning counter (+1).

The transmitPrestige [ toolType ] takes as toolType either “pot” or “point” and is used to determine which attributes are copied.

  1. Ego identifies who on the same patch they will learn from. For points, they select the person with the highest pointPrestige. If potPrestige is T, they choose the person with the highest potPrestige, but if potPrestige is F, they choose someone at random. This person is now the “teacher.”
  2. The teacher adds its list of traits to a temporary global list.
  3. If visibility affects transmission, the list is filtered through the getMostlyVisibleTraits reporter.
  4. The reporter chooses and reports one of the traits in the teacher’s list. The interaction links between the teacher and Ego updates its learning counter (+1).

The getMostlyVisibleTraits [ traitList ] reporter selects and reports mostly visible tools from the traitList given. There is some error as the visibility of a tool becomes its probability of being chosen.

  1. This first separates the traitList into two lists. One that focuses on the main types, and one that focuses on types with variants. For the main types, it simply pairs this new list with a list of 2s of the same length, as all main tools have high visibility. For variants, it pairs the list with each variants’ visibility obtained from the visibilityTable.
  2. It then merges both paired lists and samples traits from that list, using their score as the probability to be chosen. The list of chosen traits is as long as the original list and includes replacements.
  3. The reporter then removes any duplicates and reports this list.

Cultural innovation

Innovation can occur through the procedures luckyLeap (where a person creates a new main type) and createVariant (where a person creates a variant of a main type they know) (see Kolodny et al. 2015; Creanza et al. 2017 for a similar implementation of cultural dynamics in an ABM. This allows creating innovations that are dependent on other innovations, which is a key feature of cultural innovations. For example, the invention of a new snare hunting method may enable new kinds of game hunting and in turn the further invention of specialized tools for processing these novel foods. Such dependence in culture is what makes it cumulative over time.These procedures add the newly invented trait to a person’s repertoire and to the master pot and point traits tables. Main types are coded as integers, and their variants are coded as letters.

The luckyLeap [ nextTrait toolType ] procedure is called by a person (ego) in Go if a random value is below the threshold set in pLuckyLeap (Interface). The value is multiplied by 2 for pots so that pots are twice less likely to see random innovations.

The nextTrait is the number of the next trait in the list of main tools to invent, calling getNextNumber to calculate it. The getNextNumber reporter takes the list of pot or point attributes ego has and reports the number following the number ego already knows. That number is converted as a string to be used by luckyLeap. The toolType can be either “pot” or “point” and is used to determine which traits are innovated.

  1. Ego adds the nextTrait to its relevant trait list (point or pot).
  2. It then removes any duplicates.
  3. If this type is not in its relevant master table, it adds an entry for it.
  4. If toolType is “point” and the traitEnvironment table does not have it already, ego adds an entry to that table with the key = trait, and the value = the environmentZone of the patch where ego stands.

The createVariant [ toolType ] procedure is called by a person in Go if a roll die is successful. It takes one main type from the agent’s known traits and adds a new variant to it.

  1. It first populates a temporary variable with the list of traits the agent knows for the given toolType.
  2. It then chooses one trait at random from that list. It uses the getTraitNumber reporter to get the number of the trait (main type). It also identifies the last variant the agent owns of that type, using the getLastVariant reporter.
  3. The procedure then checks if ego has reached the limit of possible variants (“o”). If they have not, they identify the next variant in the alphabet. This letter becomes the newVariation.
  4. Ego merges the newVariation to the main type to create a new type with variant.
  5. If toolType is “point” and the traitEnvironment table does not have it already, ego adds an entry to that table with the key = trait, and the value = the environmentZone of the patch where ego stands.

Creating the archaeological record

The dropItem procedure removes objects from a person’s inventory based on the pointsBreakageRate or potsBreakageRate parameters, which are now set to 0.3. If the person dropping an object is hunting, the probability of discard is unchanged. If not hunting, the probability is divided by twenty to represent less likelihood of losing a point while engaged in other activities. This value was chosen to represent the probabilities of a point being discarded at a camp based on Claudine Gravel-Miguel et al. 2021. The potsBreakageRate does not vary. When dropping an item, a person removes the first item in its inventory (the oldest item made). If the person is on a camp or aggregation site, that site records the discarded object in its pot or point assemblage.

Global Extensions Description {-} These are the extensions used in the model. Note that stats is not bundled and must be downloaded and placed in the extensions folder of the Netlogo app or in the same directory as the model. sr is commented out as it requires additional setup, but can be enabled by uncommenting.

Global Variables Description {-} These are the global variables defined for the model.

Agent Types Description {-} Each agent type forms a nested hierarchy–bands contain camps which contain people.

Agent Properties Description {-} These are the variables specific to each type of agent.

Patch Variables Description {-} These variables are specific to patches.

Link Variables Description {-} These variables are specific to links.

Procedures setup Description {-} Setup the model.

setupPatches Description {-} Procedure to create patch environments (1 or 3).

setupBandsAndTerritories Description {-} Procedure to create bands and their territories. Each band is separated by a minimum distance.

setupCampsAndPeople Description {-} Procedure to set up camps and their people.

setupMasterTables Description {-} Observer procedure. Creates 3 tables, one for point traits, one for pot traits, and one for points’ environment (in which environment aare they fit?).

setupVisibility Description {-} Observer procedure. Creates visibility table.

go Description {-} This is the main procedure that calls the others. In this order: the visibility table is updated; the environment is changed (if necessary), the adjacency tables are updated, alliances are cleaned, innovations are made, people craft new object, aggregation occurs (if applicable), people do activities, prestige gets calculated, Traits are lost, traits and items are transmitted, objects are dropped, migration occurs (if applicable). Then, if the run has reached its end, metrics are outputted (if applicable).

changeActivity Parameters {-} Description needs to be a string identifying the new activity. It can be “hunting”, “learning”, “trading”, or “visiting”.

Description {-} Change the activity of the person and if recordPerson is true then records the activity.

aggregate Description {-} Observer procedure. Tells people to go to their band’s spot, which serves as the aggregation camp. All people in the same camp add 1 visiting to their counter. People can then create alliance and exchange items with people in the same camp.

disperse Description {-} Observer procedure. Called when the aggregation is over. Tells people to return to their home camp.

hunt Description {-} People procedure. People go hunting with band mates.

createAlliance Description {-} People procedure. People create a new alliance. When aggregating, they choose someone new (who is not already their ally) in the same camp. At other times, they choose someone new from another camp, favoring people from nearesst camps.

cleanUpAlliances Description {-} People procedure. People count the number of allies they have and drop one or more if they have more than the limit. Allies with strong ties have less chance to get dropped than allies with weak ties.

visit Description {-} People procedure. People go visit one of their allies.

calculatePrestige Description {-} People procedure. People take their inventory of points and calculate how many of them are fit in the current environment. This gets added to their huntResults variable. They take the average of the last 10 of those values as their prestige.

tradeItems Description {-} People procedure. Two people in the same place exchange one of their point and/or pot with one another. When they do so, they add 1 to their “trading” counter.

dropItem Parameters {-} itemType can be “point” or “pot” so it can be called individually. prob? is a boolean used to increase flexibility of the procedure (so it can be called by breaking the object as well as if they are over their inventory limit).

Description {-} People procedure. People drop one itemType with a probability set by the object’s breakageRate (if prob? is FALSE) or automatically if they are over their hand limit (prob? TRUE). Hunters are much more likely to drop a point when hunting. If it was the last point they had, they create a new one. If the drop occurs at a camp, that camp records the object’s characteristics.

returnToCamp Description {-} People procedure. Agents move back to home camp or to the aggregation location (can be called by disperse or doActivity at the end of a hunt).

updateCampInteractions Parameters {-} interactionType is a string that can be “learning”, “hunting”, “trading”, or “visiting”. agentA and agentB are two different people doing the interaction.

Description {-} Observer procedure. Update the camp interaction link every time agents interact

craftItem Parameters {-} itemType can be “point” or “pot” so it can be called individually. prob? is a boolean used to increase flexibility of the procedure (T: it happens with a certain probability, F: it happens automatically).

Description {-} People procedure. People have a certain probability to create an object (pot or point). When created, the object records the current tick as well as its type (one of the types its maker knows).

migrate Description {-} People procedure. One person chooses another camp and trade place with someone from that camp. Migration can be gendered, in which case, only people with gender 1 trade places.

recordCampPrestige Description {-} Camp procedure. Camps calculate the mean and max prestige of their people.

luckyLeap Parameters {-} objectType can be “point” or “pot”. nextTrait is the trait that will be innovated (added to the list).

Description {-} People procedure. People create a new object (point or pot) with the trait set to the given nextTrait. They then create a new object.

createVariant Parameters {-} objectType can be “point” or “pot”.

Description {-} People procedure. Add a variant to a mainType already known by ego. Then create an object.

loseTrait Parameters {-} traitsList represents the list of known trait by ego (either potTraits or pointTraits).

Return {-} Returns the trait or traits that ego remembers after forgetting one trait (and its variant if it was a main type).

Description {-} People procedure. People forget one of the traits they know, at random. If a main trait is forgotten, ego forgets its variants as well.

getNextNumber Parameters {-} objectType can be “pot” or “point”.

Return {-} Returns the max number already known + 1.

Description {-} People procedure. Called with LuckyLeap only. It finds the next number (type) to innovate.

getFitTrait Return {-} Returns one point type with a relatively high fitness, given the current environment.

Description {-} People procedure. Looks at the person’s pointTraits’s fitnesses in the current environment, and then chooses one of the traits from that list (fit traits have more chance to be selected than non-fit traits).

deconstructingTrait Parameters {-} trait is a trait known by the person calling this reporter. It can be either a point or a pot trait (not the object made with it).

Return {-} Returns the last part of the trait in its correct format (either a letter as a string or a number as a number).

Description {-} People procedure. Used to define if a trait is a main type or a type with a variant, so the person innovating onto it can add the correct variant (called in createVariant).

getTraitNumber Parameters {-} trait is a trait known by the person calling this reporter. It can be either a point or a pot trait (not the object made with it).

Return {-} Returns the number of the trait as a number (main type).

Description {-} People procedure. Gets the main type of the object (without variants).

nextLetter Parameters {-} currentLetter is the letter of the variant upon which the person will innovate.

Return {-} returns the next letter alphabetically.

Description {-} People procedure. Helps find the next letter when innovating on a trait with variant.

traitVisibility Parameters {-} trait is one trait known by ego (can be for points or pots). The trait needs to be one with variant.

Return {-} Returns the visibility assigned to the trait provided.

Description {-} People procedure. Identifies the visibility of a certain trait. Called by checkVisibility only for traits with variants.

getMostlyVisibleTraits Parameters {-} traitList can be a pointTraits list or the potTraits list.

Return {-} Returns a list of some of the most visible traits from the traitList provided.

Description {-} People procedure. Samples traits from the traitList with repeat. Visible traits have higher chances to get selected.

getVariantsOfMainType Parameters {-} trait is the trait considered, objectType can be “point” or “pot”.

Return {-} Returns a list of all the variants of the trait’s mainType already known by ego (NOT all variants created in the simulation).

Description {-} People procedure. Identifies all the variants of trait’s mainType known by ego.

getLastVariant Parameters {-} trait is the trait considered, objectType can be “point” or “pot”.

Return {-} Returns the last variant ego knows of the type of the trait provided (variant based on alphabetical order).

Description {-} People procedure. Identifies the last variants ego knows of the type of the given trait.

updateVisibility Description {-} Observer procedure. Assigns visibility of 1 (low visibility) or 2 (high visibility) to a trait or variant.

transmitCulture Description {-} People procedure. Calls the appropriate procedure to transmit traits based on the transmissionRate, learning method, and object selected.

transmitConformism Parameters {-} objectType can be “point” or “pot”

Return {-} Returns the most common trait

Description {-} People procedure. Ego identifies the most common trait known by most people at the same spot. All people who contributed add 1 to their learning counter with ego.

transmitPrestige Parameters {-} objectType can be “point” or “pot”.

Return {-} Returns the trait chosen from the person with higher prestige.

Description {-} People procedure. Ego identifies the person with highest prestige around and takes one of its traits. Ego adds 1 to its learning counter with the person with prestige.

totalTraits Parameters {-} trait can be “potTraits” or “pointTraits”.

Return {-} Returns the number of traits currently known by the population.

Description {-} Observer procedure. Called by the Trait counts plot.

printMetrics Description {-} Observer procedure. Print multiple metrics in the Command Center.

assemblageSize Parameters {-} objectType can be “pot” or “point”.

Return {-} Returns the number of objects discarded at camps and bands.

Description {-} Observer procedure. Calculates the number of objects (pot or point) discarded at camps and bands.

interactionFreq Parameters {-} activityName is a string describing the activity for which we want the interactions (hunting, learning, trading, or visiting).

Return {-} Returns the rounded sum of all interactions of the given activity.

Description {-} Observer procedure. Calculates the sum of the interactions of the given activity.

getBRBands Parameters {-} objectType can be “pot” or “point”.

Return {-} Returns the mean of all inter-band distances (Brainerd-Robinson). Scaled where 1 is most similar, and 0 is no similiarity.

Description {-} Observer procedure. Computes the Brainerd-Robinson similiarity between bands.

getBRCamps Parameters {-} objectType can be “pot” or “point”.

Return {-} Returns the mean of all inter-CAMP distances (Brainerd-Robinson). Scaled where 1 is most similar, and 0 is no similiarity.

Description {-} Observer procedure. Computes the Brainerd-Robinson similiarity between camps of the same band.

calculateBR Parameters {-} objectType can be “pot” or “point”.

Return {-} Returns a list of distances (Brainerd-Robinson) between the agent asking and its linked agents. Scaled where 1 is most similar, and 0 is no similiarity.

Description {-} Observer procedure. Computes the Brainerd-Robinson similiarity of the objectType between linked agents.

listTraits Description {-} Observer procedure. Tells bands and camps to create a list with all the traits known by their people (one list for point and one list for pot). Then collate this list into a table that has the frequency of each trait type.

getTraitFreqTable Parameters {-} allTraits is a table of all traits collected by either the camp or band. That table has the traits known by their people (key) and the count of those traits (value).

Return {-} Returns a table with the frequency of each trait given as a percentage rather than a count.

Description {-} Camp AND Band procedure. Transforms the count of each trait into a percentage of known traits (frequency).

getFreqDistList Parameters {-} traitsTableA and traitsTableB are two tables of traits with key = traits, and values = frequency of that trait known by the people of the band/camp asking.

Return {-} Returns a list of the frequency differences for all traits found in the two tables.

Description {-} Band AND Camp procedure. Calculate the difference in popularity between all items in two tables.

modularityBand Parameters {-} objectType can be “pot” or “point”.

Return {-} Returns the modularity of bands. Modularity is “defined based on the number of in-community links versus the number of between-community links” (manual).

Description {-} Observer procedure. Calculates the modularity of bands, based on the traits known. Creates links between people based on the similarity of their known traits and uses those links to calculate modularity.

modularityCamp Parameters {-} objectType can be “pot” or “point”.

Return {-} Returns the modularity of camps. Modularity is “defined based on the number of in-community links versus the number of between-community links” (manual).

Description {-} Observer procedure. Calculates the modularity of camps, based on the traits known. Creates links between people based on the similarity of their known traits and uses those links to calculate modularity.

resetTempLinksForModularity Parameters {-} objectType can be “pot” or “point”.

Description {-} Observer procedure. Calculates the Jaccard distance between traits of people. The most similar create a new set of links used to compute modularity.

networkMetric Parameters {-} objectType can be “pot” or “point”.

Return {-} Returns a value between -1 and 1. If it is close to 0, it means the network from artifacts represents relatively well the network from interactions, if it is positive, the network from interactions is more clustered than the network from artifacts, if it is negative, the network from artifacts is more clustered than the network from interactions.

Description {-} Observer reporter. Calculates the difference between the clustering coefficient of the network of people and the clustering-coefficient of the network reconstructed from similar objects.

setupJaccardDistBetweenTraits Parameters {-} objectType can be “pot” or “point”.

Description {-} Link procedure. Calculates the Jaccard distance between traits of agents linked by the link asking.

setupJaccardDistBetweenObjects Parameters {-} objectType can be “pot” or “point”.

Description {-} Link procedure. Calculates the Jaccard distance between assemblages (objects) of agents linked by the link asking.

jaccardDist Parameters {-} A and B are two lists (of traits or objects discarded) to be compared.

Return {-} Returns the Jaccard distance between the two lists.

Description {-} Observer procedure. This calculates the Jaccard distance between two lists (A and B).

jaccardDistWeighted Parameters {-} A and B are two lists (of traits or objects discarded) to be compared.

Return {-} Returns the weighted Jaccard distance between the two lists.

Description {-} Observer procedure. This calculates the weighted Jaccard distance between two lists (A and B).

PhiSTCamps Parameters {-} objectType can be “pot” or “point”.

Return {-} Returns the PhiST statistic. If it is close to 0, variation is not significantly structured around camps.

Description {-} Observer reporter. Calculate PhiST statistic between people’s inventories according to the camp they belong to.

PhiSTBands Parameters {-} objectType can be “pot” or “point”.

Return {-} Returns the PhiST statistic between bands. If it is close to 0, variation is not significantly structured around camps.

Description {-} Observer reporter. Calculate PhiST statistic between people’s inventories according to the camp they belong to.

nTraits Parameters {-} objectType can be “pot” or “point”.

Return {-} Returns the number of traits known of the given objectType. It considers repetitions (e.g., if 15 agents all know trait “1” only, the result will be 15).

Description {-} Observer procedure. Calculates the number of traits known by the whole population (with repetitions).

compileArtifactSimilarities Description {-} Observer procedure. Calculate the similarity of objects between camps, keep only the most similar and update the relevant adjacency tables with the similarity data.

compileCampInteractions Description {-} Observer procedure. Asks the links between camps to update the adjacency tables with their counters.

updateAdjacencyTables Parameters {-} matrixName is the name of the adjacency table getting updated, and addedValue is the value that will be added to the table.

Description {-} Link procedure. Updates the values in the adjacency table given (matrixName) by adding addedValue to the values already there.

fileNameBase Return {-} Returns a filename with the following structure: “day-month-year-time”.

Description {-} Observer procedure. Takes the current time and produces a name with it.

exportResults Description {-} Observer procedure. Exports variable and values to a CSV file with filename set by fileNameBase.

exportFineTunedData Description {-} Observer procedure. Writes adjacency tables data to a separate file at 1000 ticks interval if wanted.

getJaccardBand Parameters {-} objectType is either “points” or “pots”.

Description {-} Observer procedure. Calculates the Jaccard distance between the traits found in the camps belonging to each band. Essentially this is a measure where 0 is most similar and 1 is most distant for sharing the same types and variants between bands.

getJaccardCamp Parameters {-} objectType is either “points” or “pots”.

Description {-} Observer procedure. This is the same as the getJaccardBand procedure except the values are calculated between camps.

calculateJaccard Parameters {-} objectType is either “points” or “pots”.

Description {-} Observer procedure. This function calculates the Jaccard distance between two lists.

adjCompare Parameters {-} m1 and m2 are matrices.

Return {-} Returns a number between zero and one where one represents identical matrices and zero is no similarity.

Description {-} Observer procedure. This function calculates the absolute differences between two matrices. It is a method to compare networks.

adjCorr Parameters {-} m1 and m2 are matrices

Return {-} Returns a number between -1 and 1 representing the Pearson’s R correlation.

Description {-} Observer procedure. This procedure uses the sr extension which is disabled by default and this code is commented out. To enable. Uncomment the sr extension, follow the sr setup instructions, install R on your computer with the igraph and multinet packages, and uncomment this code. This procedure takes two matrices and converts them to networks. It then calculates the Pearson’s R correlation between the network edges. It is a method to compare networks.

GUI elements | Name | Type | Value | Properties | | —- | —- | —– | ———- | | nTicks | INPUTBOX | 2000.0 | entrytype: Number | | randomSeed | INPUTBOX | 9433.0 | entrytype: Number | | randomizeSeed | SWITCH | TRUE | | | saveResults | SWITCH | TRUE | | | pLoss | SLIDER | 33 | min: 0; max: 100; incr: 1 | | aggregationFreq | SLIDER | 12 | min: 6; max: 36; incr: 1 | | pNewObject | SLIDER | 10 | min: 0.05; max: 100; incr: 0.05 | | aggregationDuration | SLIDER | 1 | min: 0; max: 12; incr: 1 | | pMigration | SLIDER | 0.05 | min: 0; max: 10; incr: 0.1 | | transmissionRate | SLIDER | 5 | min: 0; max: 100; incr: 1 | | nOfBands | SLIDER | 3 | min: 3; max: 6; incr: 1 | | campPopulation | SLIDER | 15 | min: 1; max: 100; incr: 1 | | pHunting | SLIDER | 1 | min: 0; max: 100; incr: 1 | | pVisiting | SLIDER | 5 | min: 0; max: 100; incr: 1 | | huntingGroupSize | SLIDER | 5 | min: 1; max: 100; incr: 1 | | pNewAlliance | SLIDER | 0.5 | min: 0; max: 10; incr: 0.05 | | maximumAlliesN | SLIDER | 50 | min: 10; max: 150; incr: 1 | | pLuckyLeap | SLIDER | 1 | min: 0; max: 10; incr: 0.1 | | learningMethod | SLIDER | 50 | min: 0; max: 100; incr: 1 | | recordPrestige | SWITCH | TRUE | | | gradualEnvChange | SWITCH | TRUE | | | pEnvChange | SLIDER | 0.1 | min: 0; max: 20; incr: 0.1 | | uniqueTraits | SWITCH | TRUE | | | visibility | SWITCH | TRUE | | | genderedMigration | SWITCH | TRUE | | | pNewVariant | SLIDER | 6 | min: 0; max: 100; incr: 1 | | potteryPrestige | SWITCH | TRUE | | | pTrading | SLIDER | 5 | min: 0; max: 100; incr: 1 | | recordPerson | SWITCH | TRUE | | | debug | SWITCH | TRUE | | | showMetrics | SWITCH | TRUE | | | HighResolutionData | SWITCH | TRUE | |

References cited

Creanza, N., Kolodny, O., Feldman, M.W., 2017. Greater than the sum of its parts? Modelling population contact and interaction of cultural repertoires. Journal of The Royal Society Interface. 14, 20170171.

Derex, M., Mesoudi, A., 2020. Cumulative Cultural Evolution within Evolving Population Structures. Trends in Cognitive Sciences. 24, 654–667.

Hamilton, M.J., 2021. Collective computation and the emergence of hunter-gatherer small-worlds (preprint). Ecology.

Kolodny, O., Creanza, N., Feldman, M.W., 2015. Evolution in leaps: The punctuated accumulation and loss of cultural innovations. Proceedings of the National Academy of Sciences. 112, E6762–E6769.

Mesoudi, A., 2011. Cultural Evolution: How Darwinian Theory Can Explain Human Culture and Synthesize the Social Sciences. University of Chicago Press.

Migliano, A.B., Battiston, F., Viguier, S., Page, A.E., Dyble, M., Schlaepfer, R., Smith, D., Astete, L., Ngales, M., Gomez-Gardenes, J., Latora, V., Vinicius, L., 2020. Hunter-gatherer multilevel sociality accelerates cumulative cultural evolution. Science Advances. 6, eaax5913.

Powell, A., Shennan, S., Thomas, M.G., 2009. Late Pleistocene Demography and the Appearance of Modern Human Behavior. Science. 324, 1298–1301.

Premo, L.S., Tostevin, G.B., 2016. Cultural Transmission on the Taskscape: Exploring the Effects of Taskscape Visibility on Cultural Diversity. PLOS ONE. 11, e0161766.

Rorabaugh, A., 2015. Modeling Pre-European Contact Coast Salish Seasonal Social Networks and Their Impacts on Unbiased Cultural Transmission. Journal of Artificial Societies and Social Simulation. 18.

Salali, G.D., Chaudhary, N., Bouer, J., Thompson, J., Vinicius, L., Migliano, A.B., 2019. Development of social learning and play in BaYaka hunter-gatherers of Congo. Scientific Reports. 9.

Salali, G.D., Chaudhary, N., Thompson, J., Grace, O.M., van der Burgt, X.M., Dyble, M., Page, A.E., Smith, D., Lewis, J., Mace, R., Vinicius, L., Migliano, A.B., 2016. Knowledge-Sharing Networks in Hunter-Gatherers and the Evolution of Cumulative Culture. Current Biology. 26, 2516–2521.

Yasuoka, H., 2006. Long-Term Foraging Expeditions (Molongo) among the Baka Hunter-Gatherers in the Northwestern Congo Basin, with Special Reference to the “Wild Yam Question.” Human Ecology. 34, 275–296.


  1. authors listed alphabetically by first name–each contributed equally to the agent based model.↩︎