In LiYingWang/kwl-burials: Compendium of R code and data for KWL burials

knitr::opts_chunk$set(
  collapse = TRUE,
  warning = FALSE,
  message = FALSE,
  echo = FALSE,
  comment = "#>",
  dpi = 300,
  fig.path = "../figures/"
)

This supplementary document contains:

• further historical background of 17th century northeastern Taiwan in a pericolonial context
• further review of network analysis and archaeological applications relevant to this study
• further definitions of our approach to tie formation and node attributes from Kiwulan burials
• further details on our methods of modeling fitting, and
• assessment of priors, MCMC outputs and Goodness-of-fit diagnostics for our models

Kiwulan in the 17th northeastern Taiwan in a pericolonial context {-}

The involvement in long-distance exchange is often associated with changes in political-economic strategies from a corporate mode to a network mode, where ambitious individuals are able to build their local power base and personal influence through the acquisition and distribution of foreign prestige goods [@Blanton1996; @Feinman2000; @Klehm2017]. Increased social inequality associated with the use of exotic trade goods is recognized as a "network" strategy for gaining power [@Feinman2000]. Northeastern Taiwan became more involved in regional exchange with East Asia since the 14th century and long-distance trade with Southeast Asia brought by Europeans in the 17th century [@Chen2005; @Wang2007]. In 1626, the Spanish built their forts in Heping dao and Tamsui, northern Taiwan, and later were defeated and taken over by the Dutch in 1642 [@Andrade2007]. Indigenous people in northeastern Taiwan first encountered colonial power in 1632 when they were attacked by the Spanish and later in 1647 by the Dutch. The Dutch asked local settlements to pay deer hides or rice as annual tributes but negotiating a delayed payment is allowed if it is needed [@Borao2009; @Kang2016]. In the annual event held by the Dutch with participants of local Indigenous representatives [@Andrade1997], the Dutch would give those representatives symbolic items such as rattan sticks or foreign goods to assert their political authority and consolidate their relationship with Indigenous settlements. The introduction of foreign prestige goods to local Indigenous communities might have occurred before European presence became established, but was maintained and amplified after European arrival due to local competition for prestige accumulation [@Borao2009; @LiandWu2006; @LWandBM2020ornament; @Ueda2016].

load(here::here("analysis", "data", "derived_data", "burial_bergm_model.RData"))
readRDS(here::here("analysis", "data", "derived_data", "pre_bergm.rds"))
readRDS(here::here("analysis", "data", "derived_data", "post_bergm.rds"))

Review of network analysis and archaeological applications {.unnumbered}

Drawing largely on network science and graph theory, a key assumption in the construction of archaeological networks is that a phenomenon can be conceptualized as a network through abstraction and representation [@Brandes2013; @Collar2015; @Mills2017; @Peeples2019]. A social network is generally visualized as a graph consisting of a set of socially-relevant nodes/actors, connected by edges/ties representing one or more relations, such as friendship, collaborations, information flow, trade ties, or any other forms of connection of interest [@Marin2011; @Wasserman1994]. The ties can be classified into four major types, including similarities, social relations, interactions, and flows [@Borgatti2009]. In archaeology, actors can be people, groups, objects, places, or events, with ties built on similarity, proximity, or co-presence of material culture to create patterns reflecting influence, geographical distance, or affiliations in social groups [@Mills2017; @Brughmans2018; @Peeples2019]. For example, past trade can be conceptualized as a network of individual entities connected by shared similar objects or elements, such as the flow of goods, to represent their interactions [@Collar2015, pp. 4]. This concept can be applied to burial contexts where each burial is an actor linked by similar elements, such as burial goods. Burial goods, especially high value goods, can reflect social practices in broader cultural contexts to represent personal wealth or social status from which we can infer social differentiation or complexity [@Gamble2002; @Janes2013]. This enables the exploration of the structure of the past social organization through the identification of the relationships among burials.

Network analysis has been increasingly applied by archaeologists in recent years to deal with past interactions and explore the underlying mechanisms. There are two common approaches to characterize network properties at two distinct scales: node/edge level and graph level [@Peeples2019]. Node level focuses on node properties in a network, such as centrality, representing the individual influence or social prominence in a group, while graph level assesses the whole network attributes, such as density, clustering in a network, to generalize relationship patterns [@Mills2017; @Peeples2019]. By quantifying those network properties, archaeologists can answer a wide range of research questions.

Examples include exploring the political centralization in the Kofun period in Japan through the hierarchical communication network constructed by prestige goods [@Mizoguchi2013], and the investigation of long term inter-site relationships from the Epipalaeolithic to the early Neolithic in the Near East according to trade items [@Coward2013]. Regarding burials, @Sosna2013 examined spatial patterns of burials from the Early Bronze Age in Rebesovice with two hypothesized networks constructed according to cultural and chronological similarity between burials. Recently, complex network modeling has been used to evaluate networks at both node and graph level through simulations of particular processes and statistical tests of the formation of network properties [@Brughmans2013; @Brughmans2018; @Freeman2004; @Salvini2010]. Such applications include simulations and testing of food exchange modes for Ancestral Pueblos on the aggregation of households in the American Southwest [@Crabtree2015], and exploring the diffusion of fired bricks across Hellenistic Europe by comparing similarity networks of sites with random networks [@Ostborn2015]. Another example is the assessment of hunter-gatherer exchange networks structure across the Kuril Islands using bootstrap simulation based on ceramic composition [@Gjesfjeld2015].

Current approaches to network analysis used by archaeologists are mostly restricted to a single rational structure without consideration of interaction between network variables. Our use of Bayesian inference on ERGMs is the first application to archaeological data that can bring new insights to understand past social structures by characterizing network properties as a whole. ERGMs are an important family of statistical models for networks that allows direct modeling for the formation of edges, or ties, between nodes [@Robins2007]. The assumption is that possible ties in a network are random variables and dependent on actor variables or the presence or absence of other ties [@Robins2007]. In other words, networks in ERGMs are viewed as dependent variables, where network dependencies and the attributes of nodes/edges can influence the formation of a tie [@Snijders2011]. For example, nodes with similar attributes are more likely to have a relationship, such as friendship between people with the same hobby. Ties form a small structure in a network called a graph configuration, that describes the form of dependence, such as reciprocity (relationship between two actors), transitivity or clustering (relationship between two actors through a shared third actor), homophily (relationship between actors with a similar attribute), and popularity (actors have many relationships with others) [@Robins2007; @Snijders2006; @Morris2008]. Those configurations represent the structure or the property of a network and can be expressed by network statistics.

By modeling those network statistics as direct functions of ties by specifying the forms of configurations, we can generate a distribution of random networks that represent our hypothesis-based model [@Morris2008; @Robins2007]. Such a distribution consists of a large number of possible networks that enables statistical inference and comparison with an observed network [@Robins2007]. ERGMs help us understand whether an observed network shows significantly more or less of a property of interest than the random networks generated from our model assumptions.

Tie formation and node attributes from Kiwulan burials {.unnumbered}

For trade beads with substantial differences in quantities between burials, we described each burial as having one of four levels according to their distributions across all burials (Figure is in the main text). Gold-foil beads are in levels of high (>3), upper-middle (3), lower-middle (2), and low (1); carnelian beads are in levels of high (>6), upper-middle (4-6), lower-middle (3), and low (1-2); glass beads are in levels of high (>6), upper-middle (3-6), lower-middle (2), and low (1). If burial 1 and burial 2 both have high quantities of carnelian beads, then there will be a tie connecting them. The trade goods we used to define network ties are relatively durable materials, such as glass, stone, and metal, that are less affected by taphonomic factors in buried environments. In addition, each burial has a clear boundary as a single analytical unit that reduces sampling errors. Thus, we are confident that the quantities and distributions of beads in burial contexts can reflect the original differences between burials.

For burial attributes, we assigned burials to different groups according to their age, sex, burial value and the presence/absence of ritual pottery. We have three groups for the age attribute, 0-12, 12-20, and above 20; two groups for the sex attribute, male and female; and four groups for the burial value that are specified in the main paper document. Table \@ref(tab:age-sex-table) shows the number of burials, of which sex and age can be determined. Despite the exclusion of eight burials lacking clear archaeological contexts, the burials analyzed in this study consisting of 90% Kiwulan burials can effectively represent structures of burial networks.

library(flextable)
library(tidyverse)
age_sex_by_phase <-
  burial_three_period_age_tidy %>% 
  select(Age_scale, gender, Phase) %>% 
  group_by(Phase, gender, Age_scale) %>% 
  filter(!Phase %in% c("disturbed", "chi")) %>% 
  count() %>% 
  pivot_wider(names_from = c(Phase, gender), 
              values_from = n) %>% 
  rename(age = Age_scale) %>% 
  mutate(age = ifelse(age == "NA", "unknown", age)) %>% 
  mutate(age = factor(age,
                      levels = c("0-12", "12-20", "+20", "unknown"))) %>% 
  arrange(age) %>% 
  select(age, pre_male, pre_female, pre_NA, post_male, post_female, post_NA) %>% 
  mutate_all(., funs(replace_na(., 0))) 

ft_age_sex <- flextable(age_sex_by_phase) %>% 
  colformat_double(digits = 0) %>%
  set_header_labels(pre_male = "pre-European", pre_female = "pre-European",
                    pre_NA = "pre-European", post_male = "post-European", 
                    post_female ="post-European", post_NA = "post-European") %>%
  merge_at(i = 1, j = 2:4, part = "header") %>% 
  merge_at(i = 1, j = 5:7, part = "header") %>% 
  add_header_row(values = c("","male", "female", "unknown", 
                            "male", "female", "unknown"), top = FALSE) %>% 
  align(j = 2:7, align = "right", part = "header") %>% 
  align(i = 1, align = "center", part = "header") %>% 
  width(j = ~ age, width = 1) %>% 
  width(j = 2:7, width = 0.8) %>% 
  set_caption("The number of burials with identification of age and sex at Kiwulan. The unknown means no sufficient skeletal remains for determination of age or sex")
ft_age_sex

Further details on methods of modeling fitting {.unnumbered}

burn_in <- 1000
main_iter <- 40000
aux_iter <- 5000
bgof_net <-100
bgof_si <- 5000
chains <- 32

After we set our model parameters, we simulated networks in a Bayesian framework using a Markov chain Monte Carlo (MCMC) algorithm. MCMC algorithms allow estimation of posterior distributions through direct random sampling the posterior without assuming the prior comes from any specific distribution [@Hamra2013]. We obtained a posterior distribution by constructing a Markov chain that describes a sequence of moves from current state to the next state following probabilistic rules based on the approximate exchange algorithm [@Caimo2011]. This enables a random or stochastic simulation in a long run where each move does not depend on the previous move. More chains can ensure a more desirable posterior distribution that is close to the target distribution under study, or convergence. In Bayesian ERGMs, MCMC first selects a set of edges (or a set of empty pairs of actors) with equal probability, and then switches to a pair of actors at random within the chosen set [@Caimo2011]. In our case, we set the number of chains to r chains. For each chain, the number of burn-in iterations was r burn_in and the number of iterations after the burn-in was $r main_iter$. We set the number of iterations used to simulate a network at each iteration to r aux_iter.

# Figure 1 in the SI
# get the first set of diagnostic plots
x <- pre_bergm

png(filename = here::here("analysis", "figures", "008-pre-diag.png"),
    width = 10, height = 12, units = "in", res = 360)

seqq <- 4
par(mfrow = c(x$dim, 3), 
    oma   = c(0, 0, 3, 0), 
    mar   = c(4, 3, 1.5, 1))

for (i in 1:x$dim) {

  plot(density(x$Theta[, i]), 
       main = "", 
       axes = FALSE, 
       xlab = bquote(paste(theta[.(i)], " (", .(x$specs[i]), ")")),
       ylab = "", lwd = 2)
  axis(1)
  axis(2)
  coda::traceplot(x$Theta[, i], type = "l", xlab = "Iterations", ylab = "")
  coda::autocorr.plot(x$Theta[, i], auto.layout = FALSE)
  if (x$dim > 4) seqq <- seq(4, x$dim, 4)
}

dev.off()

# Figure 2 in the SI
# get the first set of diagnostic plots
y <- post_bergm

png(filename = here::here("analysis", "figures", "009-post-diag.png"),
    width = 10, height = 12, units = "in", res = 360)

seqq <- 4
par(mfrow = c(y$dim, 3), 
    oma   = c(0, 0, 3, 0), 
    mar   = c(4, 3, 1.5, 1))

for (i in 1:y$dim) {

  plot(density(y$Theta[, i]), 
       main = "", 
       axes = FALSE, 
       xlab = bquote(paste(theta[.(i)], " (", .(y$specs[i]), ")")),
       ylab = "", lwd = 2)
  axis(1)
  axis(2)
  coda::traceplot(y$Theta[, i], type = "l", xlab = "Iterations", ylab = "")
  coda::autocorr.plot(y$Theta[, i], auto.layout = FALSE)
  if (y$dim > 4) seqq <- seq(4, y$dim, 4)
}

dev.off()

Assessment of priors, MCMC outputs and Goodness-of-fit diagnostics for our models {.unnumbered}

We conducted a sensitivity analysis to understand the effect of the prior specification on our two models. To explore the sensitivity, we ran the analysis with prior values to be vague normal distribution (mean is 0 with standard deviation is 10) and a range of plausible priors (mean ranging from 0 to 5 and standard deviation ranging from 0 to 10). We found very minor differences in the posterior distributions, but the directionality was consistent for all values. The results of these are not shown here. We choose the priors that have the best autocorrelation properties. The comparison demonstrates that our results are robust to different priors that preserve the differences between our two models, especially for density, transitivity, and centralization. However, the posterior estimates in 95% confidence intervals could vary for some covariates with different priors. In general, models with our informative priors have a better convergence of the MCMC presented by the diagnostic plots.

As an informal way to diagnose model convergence, we evaluated three diagnostic visual summaries of our MCMC output, including density plots, trace plots, and autocorrelation plots for both models (Figure \@ref(fig:figure-pre-diag) and Figure \@ref(fig:figure-post-diag)) [@Hamra2013]. We can see that the diagnostic plots show a trend of stationary distributions for most of the variables. However, the degree of autocorrelation is high at 40 lag for all variables but the dyadcov variable, which indicates slow mixing and strong correlations between our variables [@Lynch2007; @Roy2020]. Then we summarize output from our two models by goodness-of-fit (GOF) diagnostics in the Bayesian framework, where the observed network is compared with the set of networks simulated from the estimated posterior distributions of the parameters of each model [@Caimo2011; @Caimo2017]. We set r bgof_net network graphs simulated from the estimated posterior distribution in ERGMs, and r aux_iter iterations used for network simulation.

# Figure 3 in the SI (code is in 007-bgof-assessment-vis.R)
knitr::include_graphics(here::here("analysis", "figures", "008-pre-diag.png"))

# Figure 4 in the SI (code is in 007-bgof-assessment-vis.R)
knitr::include_graphics(here::here("analysis", "figures", "009-post-diag.png"))

Our Bayesian GOF diagnostics summarized three distributions, including degree, minimum geodesic distance, and edgewise shared partner distributions. The diagnostics plots (Figure \@ref(fig:figure-bgof-pre-E) and Figure \@ref(fig:figure-bgof-post-E)) demonstrate that both models fit the observed networks very well for the minimum geodesic distance distribution (i.e. the number of edges between node pairs in a shortest path, @Hunter2008). The pre-European model has a better fit for the degree distribution and edgewise shared partner distribution, compared to the post-European model. Despite some observations falling outside the 95% interval for the post-European model, the fit is generally good with most observations within it. This provides a statistical approach to check how well the estimated posterior parameter distribution, based on our hypotheses, can reproduce networks with similar general structural features of the observed networks. We then compared the distribution of our observed networks, the networks before and after the arrival of Europeans, with the distribution of our hypothesized models to test an increased social inequality after the foreign contact.

knitr::include_graphics(here::here("analysis/figures/003-pre-bgof.png"))

knitr::include_graphics(here::here("analysis/figures/004-post-bgof.png"))

References {.unnumbered}

::: {#refs} :::

pagebreak

Colophon {.unnumbered}

This report was generated on r Sys.time() using the following computational environment and dependencies:

# which R packages and versions?
if ("devtools" %in% installed.packages()) devtools::session_info()

The current Git commit details are:

# what commit is this file at? 
if ("git2r" %in% installed.packages() & git2r::in_repository(path = ".")) git2r::repository(here::here())  

Word count: r wordcountaddin::word_count()

LiYingWang/kwl-burials documentation built on Aug. 29, 2021, 1:26 a.m.