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/"
)

Introduction {-}

Network analysis has been increasingly used by archaeologists to understand past relationships, interactions, and structures of observed phenomena by visualizing and analyzing relational data [@Brandes2013; @Mills2017; @Peeples2019]. Recent developments in modeling approaches enable hypothesis testing for the intercorrelations between individual elements and overall structures in networks [@Brughmans2013; @Brughmans2018; @Freeman2004; @Salvini2010]. By comparing hypothesized networks with the observed network, statistical network modeling can answer anthropological questions related to exchange [@Gjesfjeld2013; @Crabtree2015], diffusion [@Ostborn2015], or social transformation [@Mills2013] in a statistical way. Among network modeling methods, exponential random graph models (ERGMs) are stochastic models for investigating the process of network formation through dependence assumptions for relationships and simulations for network patterns [@Harris2013; @Handcock2008; @Ghafouri2020]. ERGMs are promising for evaluating the dynamic social processes behind observed archaeological networks [@Brughmans2014]. However, computational difficulties, sensitivities to uncertainties, and ambiguities in interpretation limit the practical applications of ERGMs [@Caimo2014].

In this paper, we go beyond ERGMs with a novel Bayesian approach that alleviates the computational issues and other limitations of ERGMs to enable clearer and more robust interpretations of the processes of network formation. Using burial data from an Iron Age site in northeastern Taiwan as a case study, we explore social changes by investigating the formation of material connections between burials. Social changes in Indigenous societies, when faced with colonial powers, are commonly observed in many parts of the world, especially European colonies where Indigenous economic, cultural, or socio-political aspects were substantially impacted [@Dietler2005; @Silliman2005; @Voss2005]. Recent studies demonstrate that the indirect effects of colonialism, or involvement in long-distance trade, may also have impacted Indigenous societies, this is known as a pericolonial context [@Acabado2017; @Trabert2017; @LWandBM2020ornament]. Burials are important for understanding past societies because material culture and biological records of burial behaviors can represent the social ranking of the deceased, and social relations between them [@Binford1971; @Curet1998; @Drennan2010; @Saxe1970]. Burial treatments can further indicate social complexity or inequality where foreign goods from long-distance trade were used by Indigenous people to express status [@Carter2015; @Dolfini2019].

Burial studies are usually based on characterizing each physical trace or burial variable, such as biological records, grave forms and goods, or ritual behaviors, and finally combining and comparing those individual observations to infer the organization of past societies [@Byrd1995; @Seikel2011]. Recent studies on burials using network analysis demonstrates it is a promising tool to test expectations derived from hypotheses [@Sosna2013]. Here we use network modeling to explore similarities and differences between burials by identifying underlying variable structuring burial relations. We treat burials as a complex network based on the assumption that a network is a patterned aggregation that includes individual elements (i.e. individual burials), pairwise relationships (the dyads, for example burials with similar types and amounts of grave goods), and an overall structure showing social patterns represented in the data [@Brandes2013]. Our network indicates status similarity between individuals, where a tie means two individuals sharing a similar social status or a similar social environment that higher status individuals will be given similar burial goods [eg. @Collar2015, pp. 4; eg. @Sosna2013]. In this case, network ties do not indicate exchange or transmission but rather similarities. This is because burial goods, especially high value ones, are normally cultural symbols of personal wealth or social status from which we can infer social differentiation [@Gamble2002; @Janes2013]. Our case study illustrates how a Bayesian approach to network modeling provides new insights into the formation of archaeological networks and allows testing of anthropological models to better understand dynamic social processes.

Exponential random graph models in a Bayesian framework {-}

Exponential random graph models (ERGMs) are statistical models for social networks that allow direct modeling for the formation of edges, or ties, between nodes [@Robins2007]. The application of ERGMs first appeared in archaeology with @Brughmans2014, who studied Iron age settlement patterns in Southern Spain by modeling inter-settlement visibility networks and visual control at 159 sites and at 190 sites for a later study [@Brughmans2015]. By fitting models with the observed network of archaeological data, they proposed that ERGMs are a promising approach for exploration and hypothesis testing of social processes. Similarly, @Amati2019 modeled three networks consisting of 15 sites (AD 100 to 400) in the Caribbean to explore interaction mechanisms, including proximity, inter-cultural items, and pottery types. By comparing hypothesized networks with the observed sites, they found that the presence of hub sites can be efficiently explained by multiple interdependent mechanisms instead of only one variable exclusively. Other examples, such as exploring drivers of social stability for the Valencian Bronze Age settlements in eastern Spain [@Cegielski2020], or reconstructing structurally efficient networks of Middle Bronze Age settlements in the Aegean [@Amati2018], demonstrate the usefulness of ERGMs for hypothesis testing for social processes.

However, those studies also illustrate some limitations of ERGMs, such as sensitivity to missing data and limited ability to incorporate and represent uncertainty. Also, it is difficult in ERGMs to estimate model parameters and interpret the results due to intractable likelihood normalizing constants and model degeneracy [@Caimo2014; @Jin2013]. A normalizing constant is a function of the model parameter for making probability distributions integrate to one, which becomes harder to compute with a larger set of networks [@Caimo2017bayesian]. This is also termed “doubly intractable” since both the likelihood normalizing constant and the marginal likelihood (the evidence of the posterior) are hard to derive [@Caimo2013; @Lyne2015]. Model degeneracy is another issue where probability models tend to overestimate a small number of extreme graphs by assigning too much weight, e.g., in extreme cases such as empty (all nodes unconnected) or complete graphs (all nodes connected) [@Schweinberger2011; @Caimo2014]. One solution to these limitations is to implement ERGMs in a Bayesian framework.

In this paper, we used Bayesian inference on ERGMs to test our hypotheses. ERGMs are used as the tool to provide the conceptual background for network formation and structure, allowing testing social processes through covariates. The benefit of the Bayesian part is to deal with computational issues of ERGMs by incorporating prior information about the network configurations (details in Supplementary Online Materials) to better understand dependencies of network variables and uncertainty [@Caimo2017; @Lehmann2021]. Prior information can be derived from previous data or assumptions based on the context of data. An important advantage that Bayesian modeling has over traditional ERGMs is the application of Markov chain Monte Carlo (MCMC) simulation using the approximate exchange algorithm [@Caimo2011]. With the exchange algorithm, Bayesian ERGMs avoid doubly-intractable computations by directly sampling from the not normalized part of the posterior, which alleviates the computational problems and gives better convergence results. This enables us to deal with complicated dependence patterns with ease, providing better estimations for complex social network models with heterogeneous data [@Caimo2017; @Snijders2006].

By fitting an ERGM with the approximate exchange algorithm, a Bayesian approach generates posterior probabilities that incorporate our sample data and prior information to estimate the effect of each ERGM parameter in our models [@Caimo2015; @Nemmers2019]. This allows robust and interpretable uncertainty quantification by examining the posterior distribution [@Caimo2017bayesian]. In addition, Bayesian approaches are useful to deal with missing data, which is often a problem leading to misinterpretation of networks, especially for archaeological studies. @Koskinen2010 shows that the effect of missing data can be reduced with Bayesian modeling that can accurately predict, on average, 80% of the ties when a third of data is missing.

Case Study: Exploring Social Changes in Northeastern Taiwan {-}

We use a Bayesian ERGM to study social changes in a pericolonial context at Kiwulan, an Iron Age site in northeastern Taiwan (Figure \@ref(fig:figure-KWL-map)), which was occupied from the 14th to 19th centuries. This occupation period includes the time before the European arrival, the presence of the Spanish and the Dutch in the 17th century, and finally an influx of Han Chinese in the 19th century [@Chen2007]. Our anthropological model proposes that the influence of an indirect colonial power, combined with high local values attached to imported goods, led to increased social inequality due to competition among individuals (details in Supplementary Online Materials) [@Brumfiel1994; @Clark1994]. Before the arrival of Europeans, some imported goods already entered Kiwulan because of the regional trade network with Chinese traders [@Chen2005; @Cheng2017]. We assume that the use of imported goods and their values were amplified after the European arrival. Different from Heping dao and Tamsui in northern Taiwan where European forts are located, the remote geographical setting of northeastern Taiwan resulted in less interference from European colonial powers, but more frequent trade activities due to local natural resources, such as rice and deerskins [@Chen2005; @Cheng2017]. This makes northeastern Taiwan a special area to explore the effect of European colonization in a pericolonial context that remains not well understood compared to European colonial bases [@Acabado2017; @Trabert2017].

Involvement in long-distance exchange is often associated with an increased social inequality where ambitious individuals could gain power by building personal influence through the acquisition and distribution of foreign goods, leading to social change from corporate to network societies [@Blanton1996; @Feinman2000; @Klehm2017]. Corporate-network societies represent two distinct political-economic strategies for obtaining power, viewed as a horizontal axis crosscutting traditional vertical hierarchy [@Drennan2010; @Feinman2000]. The features of network strategies include the accumulation of prestige goods, wealth differentiation, and trade monopolization through individual networks. These are contrasted with corporate strategies that stress communal ritual elements, shared power, and where wealth differentiation, if any, would be associated with corporate groups, such as age-based or kin-based, rather than individuals [@Siegel1999]. In Asia, network-based societies are frequently identified, especially for places involved in long-distance trade [@Carter2015; @Ueda2016; @Liu2006].

The observed uneven distribution of burial goods at Kiwulan may be explained as a result of unequal access to trade goods when Indigenous societies became more involved in the complex trade networks stimulated by European colonial activities in northern Taiwan. Foreign trade goods were believed to be used as prestige goods based on archaeological contexts at Kiwulan [@Hsieh2009; @Wang2011]. There are very few trade goods found in middens but more in the residential area and burials. According to Spanish historical records they are a symbol of personal influence [@LiandWu2006]. In addition, the study of spatial distribution of trade ornaments in Kiwulan residential area demonstrates a more clustered pattern after the European presence, indicative of intentional accumulation of trade goods [@LWandBM2020ornament]. The European contact and associated trade activities could be the stimulus of competition for foreign resources among aggrandizers that led to social changes.

Using the theoretical framework of corporate-network strategies, we test a hypothesis of the shift in organizational modes in relation to increased social inequality, from a more corporate to a more network society, after the European presence that can be observed in the Kiwulan burial network [@Blanton1996; @Drennan2010; @Feinman2000]. We explore whether this change can be identified through the pattern of intercorrelations of foreign prestige goods among burials. We ask: (1) Does the burial network after European arrival present a more centralized structure with fewer subgroups compared to the burial network before European contact? and if so, (2) what are the major variables affecting or forming these two networks? Can we identify the shift in variables from ritual to wealth associated with two different strategies? By answering these questions, this case study helps to expand our understanding of European colonial effects on Indigenous groups in a pericolonial context.

knitr::include_graphics(here::here("analysis", "figures", "000-KWL-map.png"))

Materials {-}

suppressWarnings(source(here::here("analysis", "code", "000-prep-data.R")))

suppressWarnings(source(here::here("analysis", "code", "001-data-tidy.R")))

We analyzed burial data collected from the excavation reports, and the original fieldwork notes for the upper component of Kiwulan (1350-1850 AD) [@Chen2007]. A total of 90 burials were unearthed from adjacent excavation squares that provide continuous stratigraphic sections suitable for temporal comparison (Figure \@ref(fig:figure-burial-map)). The preservation of human remains at Kiwulan was generally poor, so we have incomplete data for age and sex (details in Supplementary Online Materials) and lack health status data. Burials are oriented in an east-west direction on the north side of the residential area, which is indicated by post-holes and in-situ wooden posts, suggesting a well-organized spatial arrangement of houses. Previous studies report an uneven distribution of prestige goods across burials without agreement about whether this uneven distribution hints at vertical social differences. For example, @Cheng2008 interpreted the unequal distribution of glass beads, especially the gold-foil beads between burials, as evidence for hierarchy, indicating a stratified society. However, @Hsieh2012 suggested a relatively egalitarian structure based on a comparative analysis of the frequencies of all burial goods. She found that the burials with rare prestige goods were usually associated with elders, which might indicate achieved, rather than inherited, status. One important limitation of these previous studies is that they did not use analytical units suitable for comparing behaviors from different phases. Here, we adopt a new chronological framework for the burials to test if network configurations differ from the pre-European period (before European colonization) to the post-European period (during and after European colonization until the arrival of Han Chinese immigrants).

suppressWarnings(source(here::here("analysis", "code", "002-burial-location.R")))

knitr::include_graphics(here::here("analysis", "figures", "002-burial-map.png"))

# number of each phase
library(tidyverse)
pre_E <- burial_three_period_age_number %>% filter(Phase == "pre") %>% pull()
post_E <- burial_three_period_age_number %>% filter(Phase == "post") %>% pull()
chi <- burial_three_period_age_number %>% filter(Phase == "chi") %>% pull()
disturb <- burial_three_period_age_number %>% filter(Phase == "disturbed") %>% pull()
diff_pre_post <- post_E - pre_E
pre_and_post <- pre_E + post_E

# missing values
age_miss <- sum(NAs_age$n)
sex_miss <- sum(NAs_sex$n)

To compare burial networks, we assigned burials to the pre-European period (n = r pre_E), European and post-European period (n = r post_E). Our assignments are based on an established fine-grained chronology that was reexamined and cross-validated by diagnostic materials, stratigraphic data, depth, and radiocarbon ages [for more details see @LWandBM2020ornament; @LWandBM2020pot]. We excluded burials from the Chinese phase (n = r chi) due to the smaller sample size that corresponds to the population decline at Kiwulan in the Chinese period [@Chen2007]. For the two periods studied here, we evaluated differences in network size using a vertex bootstrap method. We also excluded r disturb burials that were heavily disturbed by modern construction which prevented accurate determination of their chronology. Most burials have clear stratigraphic contexts so modern or ancient disturbance is readily visible (i.e. if a burial has been cut into by later burials, very few cases of this have been documented at Kiwulan). By comparing networks from the pre-European period and the post-European period we can examine the effects of foreign contact on community relationships at Kiwulan.

Contextualizing Burial Networks {-}

We built networks where burials (nodes in the network) are linked when they have the same trade goods in common. This is based on the assumptions that social status and associated relations can be represented by sharing similar goods and burials can reflect social structures [@Binford1971; @Saxe1970; cf. @Coward2013, pp. 252]. Thus, individuals with the same prestige goods reflect similar access to trade, exchange, and gifting networks, indicative of the similarity among Kiwulan residents and social inequality. The trade goods we identified include gold-foil beads, carnelian beads, glass beads, Chinese porcelains, stonewares, gold foils, fish-shaped ornaments, and imported coins, which were also items less affected by taphonomic processes because of their durable materials. Those goods are considered as prestige goods due to their rarity in burial contexts and descriptions in historical records [@Chen2007; @Cheng2008; @Hsieh2009; @Hsieh2012; @Wang2011]. European historical accounts mentioned that those items were treated as prestige goods in Indigenous culture [@Borao2009; @LiandWu2006]. For example, a Spanish priest described a Spanish soldier exchanging carnelian beads for natural resources with local Indigenous people, because of the beads’ high value in Indigenous culture [@LiandWu2006]. In addition, Spanish visitors observed that Indigenous people possessing more imported goods were recognized as having higher status in their community [@LiandWu2006].

Methods {-}

load(here::here("analysis", "data", "derived_data", "burial_bergm_model.RData"))

Hypothesis and prediction {-}

Our social change hypothesis is that differential access to foreign prestige goods after the European presence led to increased social inequality at Kiwulan, where the social structure changed from a more corporate strategy to a more network strategy, presenting different organizational modes of societies [@Blanton1996; @Drennan2010; @Feinman2000]. If social inequality gradually increased after European arrival as we hypothesize, then we expect to observe a network structure with higher centralization and less transitivity, along with network ties determined by wealth, indicative of a network-based society. In contrast, we predict the network before European arrival will show low centralization and high transitivity, and network ties associated with ritual elements, indicative of a corporate-based society. High transitivity (or clustering) indicates multiple subgroups with shared burial goods, representing social cohesion instead of individual accumulation in an earlier time. To test this hypothesis, we use Bayesian ERGMs to model the formation of network ties and the underlying mechanisms that shaped relationships between people at Kiwulan.

suppressWarnings(source(here::here("analysis", "code", "005-network-diagrams-two-phases.R")))

ave_value <- round(mean_burial_value,1) 
high_value <- round(high_burial_value,1) 

Construction of Network Variables {-}

Trade beads are found across burials with substantial differences in quantities, so we described each burial as having one of four levels, high, upper-middle, lower-middle, and low, according to their distributions across all burials (Figure \@ref(fig:raincloud-beads); details in Supplementary Online Materials). For less frequent prestige goods, including imported ceramics, gold foils, fish-shaped ornaments, and coins, we linked two burials when they both possess each type of goods (i.e. presence or absence). Node attributes here include osteological data, such as age and sex, and cultural data, such as ritual pottery, and a burial value index. Ritual pottery was identified as locally-made ceramics placed above an individual, which is interpreted as vessels used for funeral feasting or offering [@Hsieh2009]. We calculated our burial value by summing values of all types of prestige goods in each burial context. For example, for gold-foil beads, we take the total number of burials (90) and divide by the number of burials with gold-foil beads (46), to give a value of 90/46 = 1.96. We repeat this for each type of prestige item. For each burial we sum these values of each type of prestige item to get a burial value [@Jorgensen1991].

knitr::include_graphics(here::here("analysis", "figures", "000-raincloud-beads.png"))

We then assigned the burials into three ranks according to breaks in the distribution of burial values, high (>r high_value, top 10 percent), above average (r ave_value-r high_value), and below average (<r ave_value), as an index of wealth. Since burials tend to have multiple goods in common, the network ties are weighted instead of binary (the value 1 represents a tie and the value 0 otherwise) [@Snijders2011]. For example, if two burials have both low quantity of glass beads and porcelain in common, the tie is given a value of 2. We are able to weigh ties due to durable materials of trade goods that are preserved well in each burial. Our networks are non-directed, which means a mutual relationship where the tie between any two actors is bidirectional. The networks constructed based on these principles show that the network after the European presence has more node connectivity in general with some nodes having a larger number of connections (Figure \@ref(fig:figure-network-diagrams)).

knitr::include_graphics(here::here("analysis", "figures", "005-network-diagrams.png"))

Model specification in a Bayesian framework {-}

We quantify the relations among burials and test our hypothesis of social change using the R programming language [@Rlanguage2019] with the bergm package [@Caimo2014]. Table \@ref(tab:table-ergm-terms) lists the network parameters we used for dependence assumptions that define our models [@Morris2008]. Every parameter in an ERGM has an associated algorithm for computing the probability of observing relations between two burials. Based on our hypothesis, we model a network with increased social inequality to be represented by endogenous network effects: low transitivity and high centralization. We include burial-specific attributes as covariate effects for homophily, such as age, sex, ritual activity, and the degree of wealth, to explore the importance of these variables in tie forming. For example, if age-homophily is important here, then people of the same age should have the same burial goods. We also include the physical distance between burials as an indicator of a kinship-based relations that we assume the deceased from the same family were buried nearby [@Chen2007]. Our model may reveal the emergence of social inequality via the presence of a few individuals as network centers, having more relations with others, with more wealth differentiation after the European arrival.

ergm_terms <- readxl::read_excel(here::here("analysis", "data", "raw_data", "ergm_term.xlsx")) 
ergm_terms_selected <-
  ergm_terms %>% 
  select(-"Description", -"note")
knitr::kable(ergm_terms_selected,
             caption = "Network parameters in ERGMs used for model specifications with associated archaeological interpretation for burial relations. ERGM terminology is from the R package statnet (Handcock et al., 2008; Morris et al., 2008)")

prior_mean<- data.frame(cbind(pre_prior_mean, 
                              post_prior_mean)) %>% 
  mutate(mean_same = ifelse(pre_prior_mean == post_prior_mean, pre_prior_mean, NA))

prior_sigma <- data.frame(cbind(diag(pre_prior_sigma), 
                                diag(post_prior_sigma))) %>% 
  rename(pre_prior_sigma = X1, post_prior_sigma = X2) %>% 
  mutate(sigma_same = ifelse(pre_prior_sigma == post_prior_sigma, pre_prior_sigma, NA))

priors <- cbind(prior_mean, prior_sigma)
bgof_si <-10000

Choice of prior values to evaluate our anthropological model {-}

Priors are a distinctive aspect of Bayesian methods that enable us to combine our beliefs based on prior research with newly observed data [@Gill2014]. The normal distribution N($\mu$, $\sigma^2$), where $\mu$ is the mean and $\sigma^2$ is the standard deviation, is typically used for priors when modelling network parameters because it is sensible and convenient, allowing adjustments for a wide range of prior information [@Caimo2017; @Normand1992]. The larger the mean of the prior distribution, the more positive effect of that parameter on network formation. A larger standard deviation makes the prior more vague, enabling comparison between variables for which we do not have specific assumptions. Since low network densities are commonly found in the real world [@Caimo2017], we specified the prior of the edge density parameter to low for both network models, following a normal distribution with mean at r priors$mean_same[1], and standard deviations at r priors$sigma_same[1] (i.e (N(r priors$mean_same[1], r priors$sigma_same[1]))).

To evaluate our anthropological model of increased social inequality, specifically a shift from a corporate to network society [@Drennan2010; @Feinman2000], we specified priors that are meaningful for social inequality, especially for the variables of transitivity and centralization (see our supplementary information for more details on prior selection). For the network before European contact, we set the priors to model less social differentiation and emphasize the ritual element shared in corporate groups: higher transitivity (N(r priors$pre_prior_mean[6], r priors$pre_prior_sigma[6])), lower centralization (N(r priors$pre_prior_mean[7], r priors$pre_prior_sigma[7])), and a higher covariate effect of ritual activity (N(r priors$pre_prior_mean[4], r priors$pre_prior_sigma[4])). Conversely, to model an increased social differentiation after European contact, we set the priors for the network after European arrival to lower transitivity (N(r priors$post_prior_mean[6], r priors$post_prior_sigma[6])), higher centralization (N(r priors$post_prior_mean[7], r priors$post_prior_sigma[7])), and a higher covariate effect for burial values (N(r priors$post_prior_mean[5], r priors$post_prior_sigma[5])). For covariates about biological features, such as age, sex, we specified a vague prior (N(r priors$mean_same[2], r priors$sigma_same[2])) for both models to explore their effects. Similarly, we used the same vague prior (N(r priors$mean_same[8], r priors$sigma_same[8])) for physical distance between burials, to explore whether there is kinship-based proximity, e.g. stronger correlations for shorter distances. We then examine the fit of the observed data to the two models we specified to evaluate our hypotheses.

Reproducibility and open source materials {-}

The entire R code [@Rlanguage2019] used for all the analysis and visualizations contained in this paper is included in the Supplementary Online Materials at http://doi.org/10.17605/OSF.IO/XGA6N to enable re-use of materials and improve reproducibility and transparency [@Marwick2017]. Also in this version-controlled compendium [@Marwick2018] are the raw data for all the visualizations and tests reported here and additional details of definition of our approach to tie formation, model fitting, MCMC run details, assessment of MCMC convergence and output, and goodness-of-fit diagnostics for our models. All of the figures, tables, and statistical test results presented here can be independently reproduced with the code and data in this repository. The code is released under the MIT license, the data as CC-0, and figures as CC-BY, to enable maximum re-use.

Results {-}

suppressWarnings(source(here::here("analysis", "code", "006-bergm-distribution-stats.R")))

suppressWarnings(source(here::here("analysis", "code", "007-bgof-assessment-vis.R")))

We examined the estimates from the posterior distributions to compare their differences in structure of the simulated networks (Table \@ref(tab:table-posterior-stats); Figure \@ref(fig:figure-posterior-distribution)). For nodal covariates in the pre-European model, the values of ritual element represented by pots (ritual-homophily) are mostly greater than zero. This suggests a positive effect on the formation of relations between burials, while the values of wealth rank represented by burial values (wealth-homophily) are mostly negative, indicating it has less influence on the tie formation. This demonstrates that burials with ritual pottery tend to form relations, but burials in the same wealth level tend not to. For other covariates with values mostly around zero, such as age (age-homophily) and sex (sex-homophily), we do not see any clear tendency. Most values are positive for the physical distance variable (dyadic covariate), indicating that physical proximity between burials influences similarity in burial treatment. For the endogenous network effects, transitivity presents a significant positive effect, while centralization demonstrates a negative effect. The high positive value for transitivity suggests a tendency of burials with similar burial goods to be clustered as connected communities, indicative of the presence of multiple corporate groups sharing burial goods in common. In contrast, the strong negative centralization shows there is a tendency toward decentralization that reflects most burials having a similar number of ties without any prominent burials. This might imply that individuals have equal access to trade goods in terms of the flow of goods.

# tidy data for the table of posterior stats 
posterior_esti_two_phases<-
  bergm_two_phases %>%
  pivot_longer(!phase, 
               names_to = "parameter", 
               values_to = "posterior") %>%
  mutate(parameter= case_when(parameter == "edges" ~ "Density",
                              parameter == "nodematch.age" ~ "Age-homophily",
                              parameter == "nodematch.sex" ~ "Sex-homophily",
                              parameter == "nodematch.ritual" ~ "Ritual-homophily",
                              parameter == "nodematch.value" ~ "Wealth-homophily",
                              parameter == "gwesp" ~ "Transitivity",
                              parameter == "gwdeg" ~ "Centralization",
                              parameter == "dyadcov.distance" ~ "Physical distance")) %>% 
  mutate(parameter = factor(parameter,
                            levels=c("Density", "Age-homophily", 
                                     "Sex-homophily", "Ritual-homophily",
                                     "Wealth-homophily", "Transitivity",    
                                     "Centralization", "Physical distance"))) %>% 
  group_by(parameter, phase) %>% 
  summarize(mean = round(mean(posterior), 2),
            median = round(quantile(posterior, 0.5), 2),
            `2.5%` = round(quantile(posterior, 0.025), 2),
            `97.5%` = round(quantile(posterior, 0.975), 2)) %>% 
  pivot_wider(names_from = phase, 
              values_from = c(mean, median, `2.5%`, `97.5%`)) %>% 
  select(parameter, ends_with(c("pre", "post")))

library(flextable)
# add column labels
ft <- flextable(posterior_esti_two_phases) %>% 
  colformat_double(digits = 2) %>% 
  set_header_labels(mean_pre = "pre-European", median_pre = "pre-European", 
                    `2.5%_pre` = "pre-European", `97.5%_pre` = "pre-European",
                    mean_post = "post-European", median_post = "post-European", 
                    `2.5%_post` = "post-European", `97.5%_post` = "post-European") %>%
  merge_at(i = 1, j = 2:5, part = "header") %>% 
  merge_at(i = 1, j = 6:9, part = "header") %>% 
  add_header_row(values = c("", "mean", "median", "2.5%", "97.5%",
                            "mean", "median", "2.5%", "97.5%"), top = FALSE ) %>% 
  align(j = 2:9, align = "right", part = "header") %>% 
  align(i = 1, align = "center", part = "header") %>% 
  width(j = ~ parameter, width = 1.36) %>% 
  width(j = 2:9, width = 0.67) %>% 
  set_caption("Estimated posterior means, medians, and 95% confidence intervals for each network parameter of two models. Confidence intervals that do not include zero indicate significant effects of parameters (Caimo, 2017). Posterior means close to posterior medians that indicate symmetric distributions.")
ft

knitr::include_graphics(here::here("analysis", "figures", "006-posterior-distribution.png"))

For the post-European network model, two nodal covariates, wealth rank and age, both show tendencies towards positive values in the posterior distributions. This indicates the burials in the same wealth level and age groups tend to form relations. However, there is no clear pattern for ritual and sex variables since most values are around zero. There are also no effects for the variable of physical distance. Similar to the pre-European network, transitivity demonstrates a significant positive effect, but weaker than the effect of the pre-European network. In contrast, centralization has a significantly higher positive effect than the effect of the pre-European network. This means there is a tendency toward centralization, reflecting a limited number of burials having many more ties than others. This implies that the presence of better access to trade goods and the behavior of wealth accumulation and display in burial events. In general, the post-European network model has a smaller transitivity effect and higher centralization effect than the pre-European network model. This may suggest a reduced tendency toward clustering but high tendency toward centralization after the European presence.

suppressWarnings(source(here::here("analysis", "code", "008-bootstrap-CI.R")))

One key difference between the pre-European and post-European networks is their size, with r pre_E burials compared to r post_E burials. To confirm the robustness of our comparison between these two different-sized networks, we used a vertex bootstrap technique to cross-validate the results of our Bayesian ERGMs. The vertex bootstrap is a non-parametric method that conducts resampling for all vertices (i.e. nodes) to quantify standard errors and estimate sampling variability in the network statistics of interest [@Chen2019; @Snijders1999; @Rroberts2021]. This enables the evaluation of uncertainty for networks and tests the difference between multiple networks of different sizes by examining their confidence intervals for the network population. We used the vertex bootstrap to compute endogenous network statistics, including density, centralization, and transitivity for our two networks. We explored the sample size effect by removing nodes at certain percentages (5-40%) for both networks and comparing their confidence intervals. The results demonstrate a consistent difference between the two networks at up to 35% node removals for the centralization variable as shown in Figure \@ref(fig:figure-bootstrap-CI). This shows consistency with our earlier finding of negative centralization in the pre-European period and positive centralization in the post-European period, over several conditions of node removal. Similarly consistent with our earlier result, the vertex bootstrap shows no significant difference on density and transitivity. This indicates that the network variables, including density, centralization, and transitivity, are robust to the different network sizes of our sample.

knitr::include_graphics(here::here("analysis", "figures", "008-bootstrap-CI.png"))

Our computational models reproduce networks that resemble the structural features of our observed networks according to Bayesian goodness of fit diagnostics (details in Supplementary Online Materials), showing both models fit the observed networks very well. We also compared the first three distribution moments of each observed distribution and their corresponding simulated distributions, represented by means. Figure \@ref(fig:figure-distribution-moments) shows that modeled values from the pre-European network are slightly closer to the observed values for the mean and the variance of distributions, compared to the values from the post-European network. In general, this demonstrates a slightly better fit between model and data for the pre-European network.

knitr::include_graphics(here::here("analysis", "figures", "007-distribution-moments.png"))

Discussion {-}

A striking finding in our results is the change in network properties of Kiwulan burials from a more cohesive network with multiple subgroups toward a more centralized network with concentration of connections among fewer burials after the European arrival in the 17th century. This supports our pericolonial impact model that Kiwulan changed from a more corporate-based to a more networked-based society, as indicated in the burial goods and burial attributes. Changes in mortuary practices could reflect the development of social or political complexity [eg. @Curet1998]. A corporate-based society, as we interpret the pre-European situation at Kiwulan, stresses shared power between individuals, communal rituals, and social inequality, if any, would be associated with groups [@Drennan2010; @Feinman2000; @Feinman2000political]. In contrast, a network-based society, as we interpret the post-European data from Kiwulan, presents wealth accumulation through individual networks, prestige goods manipulation, and trade monopolization [@Blanton1996]. It should be noted that the corporate-network continuum represents a dynamic process with different degrees of hierarchical complexity instead of static ideal-type stages [@Feinman2000].

Network covariates give insights into the socio-cultural factors that were associated with the acquisition and distribution of foreign prestige goods. The general pattern of posterior distributions indicates that the main factor contributing to burial networks for the pre-European network was ritual, and for the post-European network was wealth. We interpret these results as illustrating different social and economic mechanisms influencing relations between burials for the two periods. Before European arrival, ritual behaviors could have been important factors that structured the interconnections of prestige goods in the burial contexts. This might imply that ritual was a status indicator at the time before many foreign goods were introduced to Indigenous societies. After the European presence was established, status indicators shifted from ritual to wealth differences, corresponding with the characteristics of corporate/network strategies. In addition, geographical proximity of burials, implying kinship association, has some effect on sharing similar prestige items between burials in the pre-European period, but not in the post-European period. This may suggest that the use of prestige items is more individual-oriented instead of kinship-based after European arrival. Burials in the same age group tend to have similar goods after European arrival. However the sex of the buried individual did not have any effects on network formation for both networks. We should also note that their effects may be impacted by the unavailability of age (r age_miss out of r pre_and_post burials) and sex (r sex_miss out of r pre_and_post burials) values in our sample.

We argue that the social relations at Kiwulan changed because of contact with European colonial powers and associated foreign trade. This argument is supported by changes in burial network structures through time. It should be noted that manipulation of burial rituals by the living can cause a disconnect between a person’s status in life and their status represented by burial contexts [@Hodder1980; @Pearson1982]. However, this issue can be reduced by comparing the results with evidence from other archaeological contexts, such as residential areas [@Chapman2003]. Previous studies of trade ornaments from the residential area of Kiwulan suggests an uneven spatial distribution when the Europeans were active in northern Taiwan that hints at an increased social differentiation [@LWandBM2020ornament]. This is consistent with the result of a more centralized burial network after the presence of Europeans, that supports a connection between a burial and the social role of the living person. Also, ethnographic records correspond to our interpretation of burial goods as prestige goods or status items in the local Indigenous culture [@Borao2009; @LiandWu2006], so the interrelations represented by the flow of these goods observed in the burial data likely reflect social contexts of the living people at Kiwulan.

We recognize that complex local and international trade dynamics between different Indigenous groups and other foreign groups, such as the Han Chinese and Japanese, might have contributed to social changes at Kiwulan as well. However, the current changes observed from burial networks and spatial distribution of trade items in the residential area at Kiwulan demonstrate consistent changes before and after the European arrival. The frequent trade activities stimulated by the presence of Europeans led to many trade items, including Anping jars, stonewares, Chinese porcelains, and beads found at archaeological sites with Dutch occupation at different parts of Taiwan, such as Fort Zeelandia. Although shortages mentioned in Spanish documents might indicate insufficient goods were brought into Taiwan [@Borao2009], we observed a large amount of trade goods at Kiwulan. Trade goods may not have been originally owned by Europeans if we consider the pre-existent regional trade with Chinese merchants. However, the large amount of foreign goods suggests frequent and active trade activities with local places involved in a global trade network.

As a large village with continuous occupation since the 14th century in northeastern Taiwan, the archaeological evidence at Kiwulan provides a special example to explore the reaction of Indigenous people living on the periphery to European colonial presence at a settlement scale that historical records rarely offer. Was its experience with the European presence unique or typical of the region? We need more evidence from other sites in northeastern Taiwan to dive deeper into pericolonial effects, and compare evidence between sites on the periphery and colonial centers, such as Heping Dao and Tamsui. The contemporary community of descendants of Kiwulan, the Kavalan people, shows features of a corporate society, such as a matrilineal system with age group administration and ritual practices [@LiuPC2021]. This indicates that social change has been continuous in this area, from the more centralized network mode of social organization we observed in the burials after European arrival, to the current corporate model of social organization. Future work should focus on collecting data from Indigenous sites related to European contact and Chinese contact for comparison with Kiwulan.

Conclusion {-}

We presented a novel approach for studying burials to interpret social structures using ERGMs within a Bayesian framework. This case study demonstrates the methodological benefits of Bayesian inference on ERGMs to inform and enhance studies of relational data in archaeology. A Bayesian framework can better estimate archaeological data by incorporating prior information with MCMC estimations to quantify uncertainty and provide interpretable output. An important feature of Bayesian network modeling is the ability to set prior values to evaluate competing anthropological hypotheses, enabling a comparison of anthropologically-relevant models, rather than testing a hypothesis of no difference or randomness. Bayesian network modeling can be applied to a wide range of archaeological data to examine the formation of relationships using robust probabilistic inference. This enables insights into the dynamic processes of relationship formation and the underlying factors of historical trajectories of socio-cultural phenomena.

Our case study examined social changes in a European pericolonial context in northeastern Taiwan. We tested a hypothesis of changes in burial networks at Kiwulan with the evaluation of both endogenous and exogenous network effects. The results support our model that the intercorrelation between burials changed after the European colonization period in the 17th century. Before the arrival of Europeans, the burial network has a tendency of more clustered subgroups with pottery as ritual practice as the key formation mechanism. After European arrival, the network has a tendency of centralization relative to the rarity of goods. We interpret the changes in formation mechanisms of networks, observed as a change from tie formation dominated by ritual to wealth, as an increase in wealth accumulation behaviors, stimulated by the European presence and its associated long-distance trade network. This aligns with changes in horizontal hierarchy represented by corporate-network societies, that demonstrates different degrees of social inequality in relation to different strategies for achieving power [@Feinman2000]. While ethnographic records may be ambiguous on the degree of social hierarchy at Kiwulan, the archaeological record allows us to see clear signals of changes in social organization. The archaeological record is valuable here because it gives us a better understanding of the indirect European impacts in understudied pericolonial contexts of East Asia [e.g. @Acabado2017; @Trabert2017].

Acknowledgments {-}

We would like to thank the Yilan County Cultural Affairs Bureau in Taiwan for permitting access to the data used in this paper. We thank the excavation team of Kiwulan led by Yu-pei Chen for their post-excavation analysis and detailed documentation about burials that enable this study. We also thank Alberto Caimo and Steven M. Goodreau for the discussion regarding the methods at the earliest stage of analysis. This work was supported in part by the travel grant from the Department of Anthropology at the University of Washington, the predissertation summer travel grant from the Henry Luce Foundation/ACLS Program in China Studies, and the Doctoral Dissertation Fellowship from the Chiang Ching-kuo Foundation for International Scholarly Exchange (Project #DF009-A-18). We thank Ben Fitzhugh, Peter Lape, Erik Gjesfjeld, Matthew Peeples for their insightful comments on an early draft. We also acknowledge comments from anonymous reviewers that have greatly improved this manuscript.

pagebreak

References {-}

pagebreak

Colophon

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.