Appendix C — Code for Chapter 4

This chapter contains the code used for the analysis in chapter 4. The code is organized into two sections: R and R Environment. The R section contains the code for the analysis and the R Environment section contains the R environment used to run the analysis.

R

library(tidyverse)
library(magrittr)
library(Momocs)
library(ggthemes)
library(here)
library(rio)
library(igraph)
library(ggnetwork)
library(vegan)

set.seed(1010)

# library(ggimage)

# load data and remove sites with only one projectile point

sites <- import(here("data/TontoBasinNetworkSites.xlsx"), setclass = "tibble") %>%
  dplyr::ungroup() %>%
  filter(!Site %in% c("AZ V:5:130 (ASM)", "AZ U:4:009 (ASM)", "AZ U:8:318 (ASM)", "AZ V:5:121 (ASM)", "AZ U:3:128 (ASM)"))

# data frame for network nodes
nodes <- sites %>%
  select(Site, Area, Architecture) %>%
  mutate(Architecture = case_when(str_detect(Architecture, "compound, platform mound") ~ "platform mound", str_detect(Architecture, "compound, roomblock") ~ "roomblock", TRUE ~ Architecture)) %>%
  rowid_to_column("id") %>%
  dplyr::ungroup()

coords <- import(here("data/siteCoords.xlsx")) %>%
  inner_join(nodes) %>%
  arrange(id)
CeramicData <- rio::import(here("data/CeramicData.xlsx"),
  setclass = "tibble"
)

ceramicDF <- import("data/ceramicDF.xlsx")

outTonto <- readRDS(here("data/outTontoGila.Rds")) %>%
  rename(Site = ArchaMapName)

outTonto$fac <- outTonto$fac %>%
  left_join(points)

table(outTonto$pointCluster)

cTotal <- ceramicDF %>%
  group_by(Site) %>%
  summarize(Ceramics = as.integer(sum(Roosevelt, na.rm = T)))

pTotal <- table(outTonto$fac$Site) %>%
  as.data.frame() %>%
  setNames(c("Site", "Projectile Points"))

sites %<>%
  select(-any_of(c("Projectile Points", "Ceramics"))) %>%
  left_join(cTotal) %>%
  left_join(pTotal) %>%
  mutate(id = 1:n(), .before = 1)

outTonto$fac <- outTonto$fac %>%
  select(-any_of("id")) %>%
  left_join(sites %>% select(Site, id))


# rio::export(sites,here("tables/TontoBasinSites.xlsx"))

# network functions
scale_to_1 <- function(mat) {
  rng <- range(mat, na.rm = TRUE)
  (mat - rng[1]) / (rng[2] - rng[1])
}
makeGraph <- function(d, nLinks = NULL) {
  # igraph
  m <- scale_to_1(d)
  g <- graph_from_adjacency_matrix(as.matrix(d), mode = "undirected", diag = F, weighted = T)
  n_df <- tibble(id = V(g)$name) %>%
    mutate_all(as.integer) %>%
    left_join(nodes, by = "id")
  V(g)$name <- n_df$id %>% as.character()
  V(g)$Site <- n_df$Site
  V(g)$Area <- n_df$Area
  V(g)$Architecture <- n_df$Architecture
  V(g)$centrality <- eigen_centrality(g, weights = E(g)$weight)$vector

  if (!is.null(nLinks)) {
    # calculate top n (with ties) edges and convert to igraph object
    # convert to dataframe
    n <- igraph::as_data_frame(g)
    n <- bind_rows(n, n %>% dplyr::rename(from = to, to = from))
    n %<>% # sort data frame
      arrange(from, desc(weight)) %>%
      # rank by site
      group_by(from) %>%
      mutate(rank = rank(-weight)) %>%
      ungroup()
    nSub <- n %>%
      group_by(from) %>%
      slice_max(n = nLinks, order_by = desc(rank)) %>%
      ungroup()

    gRanks <- graph_from_data_frame(nSub, directed = F, vertices = n_df)
    gRanks %<>% igraph::simplify()
    V(gRanks)$centrality <- eigen_centrality(gRanks, weights = E(gRanks)$weight)$vector
    return(gRanks)
  } else {
    return(g)
  }
}

values2graphicsMod <- function(values, output = "color", palette = "Paired") {
  if (!output %in% c("color", "shape")) {
    stop("wrong parameter: output")
  }
  if (is.data.frame(values)) {
    values <- values[[1]]
  }
  types <- as.factor(values)
  num_types <- length(levels(types))
  color_map <- brewer.pal(num_types, palette)
  color_map[types]
  if ((output == "shape") & num_types > 5) {
    warning("only 5 distinct shapes available: some shapes have been repeated")
  }
  l <- vector("list")
  l[["legend.text"]] <- levels(types)
  l[["legend.pch"]] <- 21
  if (output == "color") {
    l[["legend.col"]] <- color_map[1:num_types]
    l[["color"]] <- color_map[types]
  }
  if (output == "shape") {
    l[["legend.pch"]] <- ((1:num_types) - 1) %% 5 + 21
    l[["shape"]] <- (as.numeric(types) - 1) %% 5 + 21
  }
  l
}

dists <- list()

d <- dist(coords[, c("UTMEAST", "UTMNORTH")]) %>% as.matrix()
colnames(d) <- rownames(d) <- coords$id
# convert to scaled similarity (1 highest)
d <- 1 - (d / max(d))
dists$spatial <- d

# architecture network
# convert to adjacency matrix
tmp <- sites %>%
  separate_rows(Architecture, sep = ", ") %>%
  left_join(nodes %>% select(Site, id))
m <- table(tmp$id, tmp$Architecture) %>% as.matrix()
# calculate jaccard distance
d <- vegdist(m, method = "jaccard", na.rm = T) %>% as.matrix()
# convert to jaccard similarity
d <- 1 - d
dists$architecture <- d

# ceramic network
# convert to adjacency matrix
m <- ceramicDF %>%
  select(Type, id, Roosevelt) %>%
  mutate_at(vars(Roosevelt), as.numeric) %>%
  pivot_wider(names_from = Type, values_from = Roosevelt) %>%
  column_to_rownames("id") %>%
  as.matrix()
# calculate jaccard distance
d <- vegdist(m, method = "jaccard", na.rm = T) %>% as.matrix()
# convert to jaccard similarity
d <- 1 - d
dists$ceramics <- d

# projectile point network
# convert to adjacency matrix
m <- table(outTonto$fac$id, outTonto$fac$pointCluster) %>% as.matrix()
# calculate jaccard distance
d <- vegdist(m, method = "jaccard") %>% as.matrix()
# convert to jaccard similarity
d <- 1 - d
dists$points <- d

g2List <- list()
nms <- names(dists)
for (n in nms) {
  print(n)
  g <- makeGraph(dists[[n]], nLinks = 3)
  eigen[[paste0(nLinks, "-", n)]] <- V(g)$centrality
  g2List[[n]] <- g
}

netDF <- map2_df(g2List, names(g2List), function(g, n) {
  ggnetwork(g) %>% mutate(layer = n)
})

ggplot(netDF, aes(x = x, y = y, xend = xend, yend = yend, shape = Architecture, color = Area, label = name, size = centrality)) +
  geom_edges(color = "gray", linewidth = .75) +
  geom_nodes() +
  geom_nodetext_repel(color = "black", size = 3) +
  scale_shape_manual(values = 15:18) +
  scale_color_brewer(palette = "Set2") +
  scale_size_continuous(range = c(3, 4.5)) +
  facet_wrap(~layer) +
  theme_blank() +
  theme(
    panel.border = element_rect(fill = NA),
    strip.background = element_rect(fill = "white"),
    legend.position = "bottom"
  ) +
  guides(size = "none", shape = guide_legend(ncol = 2, title = NULL), color = guide_legend(ncol = 2, title = NULL))
ggsave(here("figures/TontoNetworksMultiNet.png"), dpi = 450, width = 6.5, height = 8)
browseURL(here("figures/TontoNetworksMultiNet.png"))
paste0(nodes$id, ": ", nodes$Site) %>% paste(collapse = "; ")

eigenDF <- map2_df(g2List, names(g2List), function(g, n) {
  tibble(id = V(g)$name, centrality = V(g)$centrality, name = n)
})

meanEigenVals <- eigenDF %>%
  mutate_at(vars(id), as.integer) %>%
  left_join(nodes, by = "id") %>%
  group_by(id, name, Area, Architecture) %>%
  summarize(centrality = round(mean(centrality), 2)) %>%
  arrange(name) %>%
  pivot_wider(names_from = name, values_from = centrality) %>%
  mutate(mean = rowMeans(across(c(ceramics, points), ~.))) %>%
  mutate(mean = round(mean, 2)) %>%
  ungroup()

meanEigenVals_architecture <- meanEigenVals %>%
  select(-Area, -architecture, -spatial) %>%
  group_by(Architecture) %>%
  summarize(across(c(ceramics, points, mean), ~ mean(.x, na.rm = T))) %>%
  rename(architecture = Architecture) %>%
  mutate_if(is.numeric, round, 2)

rio::export(meanEigenVals_architecture, here("tables/TontoEigenCentrality_architecture.xlsx"))

# run all links from 3:10
# using gcor from SNA for demonstration as
# multinet package is currently unavailable as
# of the writing of this dissertation
library(sna)
allGraphs <- list()
gcor_all <- list()
for (nlink in 3:10) {
  for (n in nms) {
    print(n)
    g <- makeGraph(dists[[n]], nLinks = nlink)
    allGraphs[[paste0(nlink, "-", n)]] <- g
  }
  gcor_vars <- matrix(NA, nrow = length(dists), ncol = length(dists))
  for (i in 1:length(dists)) {
    for (j in 1:length(dists)) {
      if (i != j) {
        n1 <- allGraphs[[paste0(nlink, "-", nms[i])]]
        n2 <- allGraphs[[paste0(nlink, "-", nms[j])]]
        n1 <- igraph::as_adjacency_matrix(n1, attr = "weight") %>% as.matrix()
        n2 <- igraph::as_adjacency_matrix(n2, attr = "weight") %>% as.matrix()
        n1 <- scale_to_1(n1)
        n2 <- scale_to_1(n2)
        gcor_vars[i, j] <- gcor(n1, n2, mode = "graph")
      }
    }
  }
  colnames(gcor_vars) <- rownames(gcor_vars) <- nms
  gcor_all[[paste0(nlink, "-gcor")]] <- gcor_vars
}

gcor_all_df <- map2_df(gcor_all, names(gcor_all), function(x, y) {
  nLinks <- str_extract(y, "\\d+")[1]
  x %>%
    as.data.frame() %>%
    rownames_to_column("network") %>%
    pivot_longer(-network, names_to = "network2", values_to = "gcor") %>%
    mutate(network = str_remove_all(network, "[0-9]|\\-")) %>%
    mutate(network2 = str_remove_all(network2, "[0-9]|\\-")) %>%
    mutate(nLinks = nLinks)
}) %>%
  filter(!is.na(gcor)) %>%
  rowwise() %>%
  mutate(network = paste(sort(c(network, network2)), collapse = " -> ")) %>%
  ungroup()

ggplot(gcor_all_df, aes(x = gcor, y = network, color = network)) +
  geom_boxplot(fill = "#D3D3D3", alpha = 0.66) +
  scale_color_viridis_d() +
  theme_gdocs() +
  theme(
    legend.position = "bottom",
    legend.title = element_blank(),
    panel.spacing.x = unit(0.5, "cm")
  ) +
  labs(x = "Pearson Correlation", y = "") +
  guides(color = "none")


eigenDF_all <- map2_df(allGraphs, names(allGraphs), function(g, n) {
  tibble(id = V(g)$name, centrality = V(g)$centrality, name = n)
})

meanEigenVals_all <- eigenDF_all %>%
  mutate_at(vars(id), as.integer) %>%
  left_join(nodes, by = "id") %>%
  group_by(id, name, Area, Architecture) %>%
  summarize(centrality = round(mean(centrality), 2)) %>%
  arrange(name) %>%
  mutate(name = str_remove_all(name, "[0-9]|\\-")) %>%
  ungroup() %>%
  mutate(id = factor(id)) %>%
  filter(name %in% c("ceramics", "points"))

ggplot(meanEigenVals_all, aes(x = id, y = centrality, color = Architecture)) +
  geom_boxplot(fill = "#D3D3D3", alpha = .5, linewidth = 0.6) +
  facet_wrap(~name) +
  scale_color_viridis_d() +
  coord_flip() +
  theme_gdocs() +
  theme(
    strip.background = element_rect(fill = "transparent"),
    strip.text = element_text(size = 20),
    legend.position = "bottom",
    legend.direction = "horizontal",
    legend.title = element_blank(),
    panel.spacing.x = unit(0.5, "cm")
  ) +
  labs(x = "Site", y = "Eigenvector Centrality")

# ggsave("figures/TontoEigenCentrality_boxplot.png", dpi = 450, width = 6.5, height = 4.5)
# browseURL(here("figures/TontoEigenCentrality_boxplot.png"))


#### export session info  ####
# Get the R version
r_version <- paste0("R version=", 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_environment3.qmd")

R Environment

R version=44.3. gt=0.11.1 flextable=0.9.7 NetworkDistance=0.3.4 vegan=2.6-10 lattice=0.22-6 permute=0.9-7 ggnetwork=0.5.13 igraph=2.1.4 rio=1.2.3 here=1.0.1 ggthemes=5.1.0 huxtable=5.5.7 kableExtra=1.4.0 Momocs=1.4.1 magrittr=2.0.3 lubridate=1.9.4 forcats=1.0.0 stringr=1.5.1 dplyr=1.1.4 purrr=1.0.4 readr=2.1.5 tidyr=1.3.1 tibble=3.2.1 ggplot2=3.5.1 tidyverse=2.0.0