Simulating Geographical Variation in Material Culture: Were Early Modern Humans in Europe Ethnically Structured?

Abstract

A high degree of structuring is seen in the spatial distribution of symbolic artefact types associated with the Aurignacian culture in Upper Palaeolithic Europe, particularly the degree of sharing of ornament types across archaeological sites. Multivariate analyses of these distributions have been interpreted as indicating ethno-linguistic differentiation (Vanhaeren and d’Errico 2006), although simpler explanations such as isolation-by-distance have not been formally discounted. In this study we have developed a spatiotemporally explicit cultural transmission simulation model that generates expectations of a range of spatial statistics describing the distribution of shared ornament types. We compare these simulated spatial statistics to those observed from archaeological data for Aurignacian Europe—using Approximate Bayesian Computation—in order to test and compare a range of hypotheses concerning group interaction dynamics for the period. Among the set of hypotheses examined, we include ones where material culture does or does not drive group interaction dynamics.

Appendices

Appendices

8.1.1 Appendix 1: Bayesian Inference and Approximate Bayesian Computation (ABC)

Bayesian inference is a branch of statistics that uses observations of particular datasets to infer the probability that a proposed hypothesis, or a parameter of that hypothesis, is true. To do this, various models with set numbers of parameters are proposed, and the posterior probability distributions of these parameters are inferred using information from prior probability distributions of the parameters and information provided by the observed data, through implementing Bayes theorem. Bayes theorem states that, given parameter (or set of parameters) θ and observed dataset D, the posterior distribution of θ, denoted P(θ|D), is proportional to the product of the probability of observing dataset D given model with parameter θ, denoted P(D|θ), and the likelihood of θ, denoted π(θ), which is the distribution of θ prior to any observations being made. Mathematically, this can be written as:

$$ P\left(\theta \Big|D\right)\propto P\left(D\Big|\theta \right)\cdot \pi \left(\theta \right). $$

(8.14)

Since the explicit form of the likelihood P(D|θ) is difficult to compute in many complex problems, a family of Bayesian methods, referred to as Approximate Bayesian Computation (ABC), which do not require the likelihood function to be theoretically specified, are used (Tavare et al. 1997; Fu and Li 1997; Beaumont et al. 2002; Bertorelle et al. 2010).

In ABC techniques, a large number of datasets are simulated under a model assuming different, randomly chosen, parameter values from within prior ranges, and appropriate summary statistics are used to measure the extent to which the simulated datasets emulate the observed data. Parameter values under which the model generates datasets closest to the observed data are retained in the posterior probability distributions of the parameters.

To be able to compare the observed and simulated datasets, robust statistics that sufficiently describe the full properties of the data are used. These are called summary statistics and those developed for the current framework are discussed in detail in Appendix 3: Summary Statistics. By comparing summary statistics calculated for each simulated dataset to those for the observed data, we are able to accept to the posterior those simulations with summary statistics sufficiently close to the summary statistics for the observed dataset, referred to as the target summary statistics. The similarity δ between observed data, S, and simulated data, S’, is calculated as the sum of normalised Euclidean distances of individual summary statistics:

$$ \delta \left(S,{S}_j^{\prime}\right)=\sqrt{{\sum}_{i=1}^n\frac{{\left({s}_i - {s}_{ij}\prime \right)}^2}{{\sigma \left({s}_i\prime \right)}^2},} $$

(8.15)

where s and s’ are values of each of the summary statistics for the observed and simulated datasets, respectively, subscript i denotes the ith of n statistics, subscript j denotes the jth of N simulations and σ(s _i ’) is the standard deviation of the ith statistics over all N simulations. In performing the data analysis, we regard the ε quantile of the distribution of distances between the observed and simulated data, δ(S, S _j ’), as the best simulations—those generating data most similar to the observed data.

8.1.2 Appendix 2: Approximate Bayesian Computation (ABC) Algorithm

Let M denote the chosen model and the set of parameters of M be θ = (θ ₁ , … , θ _m ). Let S = (s ₁ , … , s _n ) and S’ = (s ₁ ’, … , s _n ’) denote the values of the summary statistics for the observed and simulated datasets, respectively. Values S = (s ₁ , … , s _n ) are referred to as the target values for each of the summary statistics. The ABC algorithm is applied as follows:

• 1. Define a set of summary statistics that capture relevant information contained in the observed dataset.

• 2. Compute summary statistics values S = (s ₁ , … , s _n ) for the observed dataset—these are the target values.

• 3. Sample parameters θ* = (θ ₁ *, … , θ _m *) from an appropriate prior distribution.

• 4. Simulate data by using parameter θ* set with model M.

• 5. Compute summary statistics values S’ = (s ₁ ’, … , s _n ’) for the simulated data.

• 6. Compute δ(S, S’), where δ is an appropriately chosen distance measure.

• 7. For a chosen tolerance ε, retain parameter set θ* in the posterior distribution of θ if δ(S, S’) < ε.

• 8. Repeat steps 1–7 until the desired number of parameter values have been accepted to the posterior distribution.

In order for ABC methods to be effective, appropriate summary statistics that sufficiently describe the observed dataset need to be developed and appropriate choices for the distance measure, δ, and tolerance, ε, must be made.

8.1.3 Appendix 3: Summary Statistics

As explained previously, to be able to compare simulated and observed datasets using ABC methods, summary statistics that capture the information contained in the observed data must be developed. These should be robust statistics and should describe sufficiently the full properties of the observed dataset considered. For the current dataset, these are:

• shared information between bead types and sites, respectively

• mutual dependence between bead types and sites, respectively

• diversity in the number of occurrences of different bead types

• cultural diversity of sites as represented by the variation in the number of distinct bead types recovered from each sites

• spatial distribution of sites

For each of these statistics, we consider the values of the mean and variance in the data analysis.

8.1.3.1 Shared Information (SI)

Shared information, denoted SI, is a statistic that measures the extent of similarity between two variables. For measuring the shared information between bead types, SI is defined to be:

$$ SI\left({t}_i,{t}_j\right)=\frac{f_i{f}_j}{{\overline{f}}^2} \log \frac{r\left({t}_i\right)+r\left({t}_j\right)}{r\left({t}_i,{t}_j\right)}, $$

(8.16)

where r(t _i ) and r(t _j ) denote the ratio of the number of occurrences of bead types i and j to the total number of sites, r(t _i , t _j ) is the ratio of the number of concurrent occurrence of bead types i and j to the total number of sites, f _i and f _j represent the number of sites in which bead types i and j occur, respectively, and $\overline{f}$ is the average number of times any bead type occurs over all sites. In this case, SI measures the similarity between pairwise bead types in terms of which sites the are present in. When two bead types never occur in the same site,

$$ r\left({t}_i\right)+r\left({t}_j\right)=r\left({t}_i,{t}_j\right),\;\mathrm{and} $$

(8.17)

$$ SI\left({t}_i,{t}_j\right)=0. $$

(8.18)

A similar equation can be used to measure the shared information between sites:

$$ SI\left({s}_i,{s}_j\right)=\frac{g_i{g}_j}{{\overline{g}}^2} \log \frac{r\left({s}_i\right)+r\left({s}_j\right)}{r\left({s}_i,{s}_j\right)}, $$

(8.19)

where r(s _i ) and r(s _j ) denote the ratio of the number of sites in which bead types i and j occur to the total number of bead types, r(s _i , s _j ) is the ratio of the number of sites that share bead types i and j to the total number of bead types, g _i and g _j represent the total number of bead types present in sites i and j, respectively, and $\overline{g}$ is the average number of bead types occurring per site. In this case, SI measures the extent of similarity between pairwise sites in terms of bead types present in those sites. Similarly to above, if two sites have no bead types in common,

$$ r\left({s}_i\right)+r\left({s}_j\right)=r\left({s}_i,{s}_j\right),\;\mathrm{and} $$

(8.20)

$$ SI\left({s}_i,{s}_j\right)=0. $$

(8.21)

8.1.3.2 Mutual Information (MI)

The mutual information, MI, between two random variables X and Y is a measure of the mutual dependence between them. It is defined as:

$$ MI\left(X;Y\right)=\sum_{y\in Y}\sum_{x\in X}p\left(x,y\right) \log \frac{p\left(x,y\right)}{p_1(x)+{p}_2(y)}, $$

(8.22)

where p(x,y) denotes the joint probability of x and y (the probability of x and y occurring together), and p ₁ (x) and p ₂ (y) denote the marginal probabilities of x and y respectively (the probabilities of the specified values of x and y occurring).

For the observed dataset in this study, setting $X={t}_i$ and $Y={t}_j$, where t _i and t _j correspond to the number of occurrences of bead type i and j in all sites respectively, allows the mutual information between all pairs of bead types to be computed. Analogously, setting $X={s}_i$ and $Y={s}_j$, where s _i and s _j correspond to the total number of bead types present in sites i and j respectively, allows the mutual information between all pairs of sites to be computed.

In contrast to the SI statistic, which only examines the common presences between sites or bead types, the MI statistic examines both the common presences and common absences. It therefore represents the dependence between the pairwise vectors in question.

8.1.3.3 Mean Absolute Deviation (MAD)

The observed dataset shows large fluctuations both in the number of bead types recovered at individual sites, and the number of times each particular bead type occurs, as shown in Fig. 8.3. Assuming that this is not the result of archaeological bias, these differences could be attributed to cultural wealth at sites, and the preference for particular bead types, respectively. To quantify this, the median absolute deviation statistic, MAD, is used. It is a measure of the variability of a random sample, and is defined to be:

Fig. 8.3

Histograms of the number of occurrences of bead types (left) and number of distinct bead types recovered from individual sites (right) for the observed data

$$ MAD= median\left(\left|{X}_i-{median}_j\left({X}_j\right)\right|\right). $$

(8.23)

Letting ${X}_i=T=\frac{f_i}{\overline{f}}$, where f _i represents the number of sites in which bead type i occurs and $\overline{f}$ is the average number of times any bead type occurs over all sites, the MAD statistic is a measure the variability in the number of occurrences of bead types. This can be thought of as a measure of variability in the popularity of, or preference for, bead types.

Letting ${X}_i=S=\frac{g_i}{\overline{g}}$, where g _i represents the total number of bead types present in site i and $\overline{g}$ is the average number of bead types occurring per site, the MAD statistic measures the variability in the number of beady types recovered. This can be thought of as a measure of variability in the cultural wealth recovered from sites.

8.1.3.4 Spatial Distribution of Sites (DR)

The extent to which sites share bead types may be a function of the distance between those sites. It is logical to expect that sites which are located near to each other share bead types more frequently than those which are far apart. The spatial distribution of sites can be explored by considering the average distance between sites sharing bead type i, ${\overline{d}}_i$, in relation to the average distance between all sites, $\overline{d}$, as follows:

$$ {DR}_i=\frac{{\overline{d}}_i}{\overline{d}}. $$

(8.24)

DR therefore quantifies the spatial distribution of sites in terms of the shared bead types between them. Figure 8.4 shows density plots for the original observed dataset (top) and a random permutation of the same (bottom). The obvious shift to the right in the density plot of the permuted dataset implies that the distance between sites sharing a particular bead type is on average larger if bead types are randomly assigned to sites. For the original observed dataset this implies that sites which are located closer to one another on average share bead types more frequently with each other than with sites that are further away, as expected.

Fig. 8.4

Density plots of the DR statistic for the original observed data (top) and a random permutation of the same (bottom)

8.1.4 Appendix 4: Bayes Factors for Model Comparison

Another useful feature of the ABC approach is the ability to formally compare the performance of different models using Bayes Factors (Kass and Raftery 1995). A Bayes Factor is a summary of the evidence provided by the data in favour of one model over another. Given models M ₀ and M ₁ , not necessarily with the same number of parameters, Bayes Factor B is given by:

$$ B=\frac{P\left({M}_1\Big|D\right)}{P\left({M}_0\Big|D\right)}=\frac{P\left(D\Big|{M}_1\right)\pi \left({M}_1\right)\;}{P\left(D\Big|{M}_0\right)\pi \left({M}_0\right)}, $$

(8.25)

where π(M _i ) is the prior probability of model M _i , P(D| M _i ) is the probability of data D given model M _i and P(M _i |D) is the posterior probability of the model, defined as:

$$ P\left({M}_i|D\right)=\frac{P\left(D\Big|{M}_i\right)\pi \left({M}_i\right)\;}{P(D)}, $$

(8.26)

where P(D) is the unconditional marginal likelihood of the data.

This form of model comparison is independent of the parameters for each model, and instead calculates the probability of the model considering all possible parameter values. This method automatically and correctly penalises model complexity; for models with a large number of parameters there is a larger parameter space to explore and so it is more difficult to find those parameter sets that generate data similar to the observed data. Therefore, models with more parameters are penalised for the increased complexity compared to simpler models, resulting in a comparison weighted by model complexity.