## Abstract

We present a new statistical approach to analysing an extremely common archaeological data type—potsherds—that infers the structure of cultural relationships across a set of excavation units (EUs). This method, applied to data from a set of complex, culturally heterogeneous sites around the Mandara mountains in the Lake Chad Basin, helps elucidate cultural succession through the Neolithic and Iron Age. We show how the approach can be integrated with radiocarbon dates to provide detailed portraits of cultural dynamics and deposition patterns within single EUs. In this context, the analysis supports ancient cultural segregation analogous to historical ethnolinguistic patterning in the region. We conclude with a discussion of the many possible model extensions using other archaeological data types.

## 1. Introduction

A central goal of archaeology is to infer the variety of cultures that were present at a given location and their relative succession in time. These estimates preferably come with corresponding ages and take account of the fact that such cultures are complex analytical entities. Traditionally, analysis towards this aim has involved inspection of diagnostic features, with fairly straightforward applications of classical statistics, such as experimental design and tests for significance [1,2]. Stimulated by analysis techniques for contemporary genetic and linguistic data, the discipline has begun to integrate the power of computational statistics [3–6]. In this work, we continue this trend by showing how statistical mixture models can be applied to one of the most common archaeological data—pottery fragments, or potsherds—to infer ancient cultural dynamics.

We develop this method in the context of a dataset comprised of potsherds derived from 11 Neolithic and Iron Age sites in the southern Lake Chad Basin of Central Africa, and classified using standard anglophone descriptive nomenclature [7]. Ethnolinguistic diversity among the populations in this area (figure 1) is among the highest in Africa and comparable to anthropologically well-known regions like highland New Guinea. This diversity is particularly great in and around the Mandara Mountains on the Cameroon–Nigeria border, occupied at high population densities by Chadic-speaking *montagnard* farmers. The historical roots of that diversity are complex: after the retreat of Lake Mega-Chad populations speaking ancestral Chadic and Nilo-Saharan languages reoccupied the region in the Mid-Holocene [5,9,10] with associated sedentism, population growth, and cultural diversification on the plains to the south between 4000 and 1200 before present (BP). Progressive occupation of the nearby Mandara Mountains by plains groups occurred as political conflict and slave raiding became regionally significant [11], with the deployment of linguistic diversity and multilingualism as sociolinguistic strategies in a challenging mountain environment [12].

Between 1984 and 2003, the Mandara Archaeological Project investigated the material correlates of this modern ethnolinguistic diversity [13–15]. As throughout much of Africa, ceramics in the Mandara Mountains are a vital material domain for the expression of identities: domestic pottery, elaborated both in morphology and in decoration, is ubiquitous and central to social interactions across societies [16,17]. In the modern Mandara context, there are strong connections between ceramic variability and the ways its production and use are associated with identity (figure 2). Simple concordance may be confounded by the sharing of ceramic traditions among groups [18], possibly due to women frequently marrying across ethnolinguistic boundaries and carrying their habits of pottery production with them [19–21]. Nonetheless, we expect that ceramic variability still largely reflects cultural traditions and long-standing social interactions both at single sites and across larger geographical scales.

Archaeological research indicates that modern Mandara populations appear to be related historically and culturally to Iron Age and Neolithic communities living on the plains adjacent to the mountains [11]. These sites date to between approximately 3500 BP and the historical period. It is likely that sites from at least the last 2 000 years are associated with Chadic-speaking farming populations. As in many equivalent Neolithic and Iron Age West African sites, pottery was recovered in very large quantities, with 312 000 potsherds comprising well over 99% of artefact inventories. As a result of the conditions of preservation of this material, very few whole pots were excavated and most potsherds are quite small. Consequently, morphological data are often not available and the primary attribute available for analysis is external decoration. As the modern context suggests a strong coupling of ceramic variability and identity, we assume that this is also true within the archaeological data. The statistical approach here then seeks to identify patterns of similar decoration distribution to infer ancient group identities.

A long-standing archaeological dilemma is how to address the incomplete mapping of physical artefacts on to cultural associations, as important shifts may not be detected by a single data type [22]. We consequently employ the term ‘cultural period’ (CP) to denote a consistent distribution of potsherd decoration types. This means that CP distributions may shift due to migration, marriage patterns, trade, other social interactions, or technological diffusion, rather than necessarily signifying the emergence of new cultures.

The core of the analysis is the Dirichlet-multinomial (DM) mixture model, similar to those used for understanding microbial ecologies [23]. In this context, the DM distribution models the frequency of decoration types within samples according to their corresponding CP. Instead of using a fixed number of CPs for the mixture, we use a Dirichlet process mixture model (DPM), a standard Bayesian non-parametric approach to mixture modelling, to integrate over all possible values [24–26]. These models are widely used in statistics, computer science, and biology in contexts where both the underlying number of model components and their parameters are unknown [27–31]. The Bayesian DPM statistical framework provides a means to infer the CPs and their underlying parameters from the data, giving a generative model of the cultural dynamics underlying potsherd production across sites and depths.

Finally, we note that the process of achaeological excavation necessitates a choice about how to best group the data for analysis by the model. Most of the field sites are subdivided into excavation units (EUs) that are in turn broken down into metre-square-sized recording units (RUs) for material recovery. At each 10 cm depth increment within a square, the removed archaeological material is collected for sorting, cleaning, classification, and further analysis; consequently, this is the most particularized form that the decoration data can take. On the one hand, data aggregated at the RU level often contained too few potsherds to provide much inferrential power, while aggregations at the site level often masks interesting variation. The analysis presented here is based on the EU-level aggregations, as this preserved the necessary variation while providing sufficient counts for the algorithm to work reliably. The data for each depth we call unit-levels (figure 2). The output of the algorithm we refer to as a culture painting.

## 2. Data and methods

### (a) Data and notation

The complete dataset is made up of 312 070 archaeological potsherd observations, collected over 31 years of field collection at 80 dig sites in the Mandara region to the south of Lake Chad on the Cameroon–Nigerian border [13]. Following standard classification procedures, these potsherds were categorized according to external decoration type, potsherd type, and interior decoration type, along with estimates of rim diameter and vessel type, where possible [32]. The statistical method presented here uses only external decoration data for three reasons: these were available for the vast majority of records; they exhibited significant variation across and within sites; and they provided significant information to previous researchers.

The dataset contains 86 external decoration types, each possessing between one and 72 100 records, with 361 associated subtypes. We filtered the dataset to contain only decoration types with at least 2 000 observations to ensure sufficient frequency for statistical inference. We also removed sites with fewer than 100 observations, leading to 11 final sites. The final dataset includes 239 629 potsherds (76.8% of the total). The final set of sites are listed according to the site designation system used for the wider archaeological project: sites 602, 618, 635, 636, 642, and 675 for those in Cameroon and 744, 755, and 756 for those in Nigeria. Most sites (7/11) contain multiple EUs, that in turn often include multiple RUs. For instance, site 642 contains 14 EUs, each of which included multiple RUs. Within each EU, stratigraphic levels correspond to depths of sediment removed for analysis, measured in 10 cm intervals beginning from the ground surface, with the total number of levels ranging from 0 to 33. Correspondingly, each potsherd is then indexed by its unit *i* = 1, … , 36, level *j* = 1, … , 34, and decoration type *d* = 1, … , 13. We use *s _{ij}* to refer to (

*s*

_{ij}_{,1}, … ,

*s*

_{ij}_{,13}), the collection of external decoration counts for unit-level (

*i*,

*j*).

### (b) Model

We take a Bayesian approach to inferring the number of CPs, their associated parameters, and their assignment to unit-levels from the data. The Bayesian perspective involves building a statistical model with two components, a likelihood and prior distribution, that specify the probabilistic relationship between the data and the model parameters and the prior beliefs about the model parameters, respectively [33]. The posterior distribution—the distribution of the model parameters given the prior distribution, the likelihood model, and the observed data—can then be derived using Bayes Theorem, as described in the Inference subsection.

For the likelihood, we adopt a mixture model formulation that assumes that the data within each unit-level come from an unobserved CP that governs the production of potsherd decorations. We assume CP *k* to have the form of a DM distribution parametrized by *A _{k}* = (

*α*

_{k}_{,1}, … ,

*α*

_{k}_{,13}), where

*α*

_{i}_{,k}is the concentration parameter for the

*i*th decoration type. A latent variable

*c*associates each unit-level (

_{ij}*i, j*) to a CP. Conditional upon

*c*=

_{ij}*k*, the likelihood at the unit-level is then: 2.1where and This likelihood allows us to account for greater dispersion in the data than would be possible with a strict multinomial distribution. The unit-levels are conditionally independent given their CP assignment and the CP parameters so it remains only to specify their prior distributions to complete the model.

The prior distribution requires we specify our existing knowledge for the number of CPs, their unit-level assignment, and their parameters (the *A _{k}*'s). Because the total number of CPs is not known

*a priori*, we employ a DPM, a non-parametric approach that provides a vague prior distribution on all possible partitions of the data into distinct CPs. This approach has two distinct advantages over a finite mixture model: it provides a full posterior distribution on the number of CPs and it shows where the model requires additional CPs to explain distributions within specific unit-levels. While both of these features could be achieved using a finite approach (for instance, using a reversible jump methodology), the non-parametric approach provides these features in a straightforward way.

Following [34], the DPM with DM components can be formulated as:
where *G*_{0} is the base measure, *γ* > 0 is a concentration parameter, denotes the likelihood in equation (2.1), and DP is a Dirichlet process. The base measure *G*_{0} is the product of two independent distributions, with the first component an exponential distribution with mean one and the second a uniform Dirichlet distribution, that is chosen for partial conjugacy with the DM likelihood [35].

### (c) Inference

We employ standard Markov chain Monte Carlo (MCMC) approaches to sample from the posterior distribution of the DPM to infer the number of CPs, the parameters *A _{k}* associated with each CP, and the assignment of unit-levels to specific CPs [34]. The implementation uses established algorithms for inference under DPMs, together with an efficient Gibbs sampling routine specific to the DM distribution [35]. Because the assignment of unit-levels to CPs is not unique up to labelling, we employ a modified version of a Kullback–Leibler minimization method to re-label the posterior samples to align CPs for each iteration with the maximum-likelihood iteration [36]. Ten independent runs of the algorithm with one million iterations provide nearly identical results, both in terms of number of cultures and inferred DM distributions. The run presented in the Results section has an effective sample size of 265.31, indicating sufficient sampling for reliable inference. A complete outline of the computational approach is outlined in the electronic supplementary material. Along with the cleaned dataset, an implementation of the model in the R computing environment is available at the Digital Archaeological Record under a Creative Commons License. The document specifics are available in the electronic supplementary material [37].

## 3. Results

Figure 3 shows the culture painting for the Mandara dataset. This presentation uses 14 primary CPs identified via an incidence plot (electronic supplementary material, figures S1 and S2), that we enumerate C1–C15, with CP 15 representing an aggregation of all clusters outside of the primary 14. Colouring indicates the fraction of posterior samples assigned to a CP over the posterior sample. Electronic supplementary material, figure S3 shows the inferred potsherd distribution for each culture and a hierarchical clustering of their relatedness, with darker shading equating to more data. Electronic supplementary material, figure S4 shows figure 3 without shading.

The culture painting shows signficant variability in CP assignment across sites, EUs and depth, as would be anticipated by previous regional archaeological and linguistic analysis. This heterogeneity in CPs across sites contrasts with strong consistency in CP assignment within the same site (e.g. sites 675 and 756) that conveys settlement continuity, while long runs of identical CPs across successive depths with occassional shifts (for instance, site 602 EU 2 and site 642 EU 2) indicate discernible shifts in production patterns. The observation of nearly all CPs across multiple sites (only CPs 2 and 3 are seen solely within a single site) and successive depth indicate the model possesses sufficient power to infer similar cultural patterns across space and time.

### (a) Culture painting preserves geographical separation

The physical separation of the western sites (602, 618, 744, 755, and 756) compared to the more easterly sites (631, 635, 636, 675, and 678) is largely preserved in the CP assignments, with the western locations dominated by CPs 1, 5, 8, 9, 10, and 13 with smaller contributions from CPs 2, 3, and 14, while the main sites largely exhibit CPs 1, 3, 4, 6, 7, 8, 10, 11, 12, and 14 with smaller amounts of CP 9. The overlapping CPs (8 and 10) between these two collections of sites largely appear at greater depths within the easterly sites. This segregation supports a long-standing geographical/cultural discontinuity between western and eastern parts of the research areas, over the last 2 000 years at least.

### (b) Distinction of Neolithic sites

Through radiocarbon dates (RCDs) and ceramic affiliations, sites 618 and 756 were previously established to be of the Neolithic period and are associated predominantly with CP 14, with site 756 nearly exclusively so [38]. The potsherd distribution for CP 14 is one of the more distinct of those inferred, as would be expected for earlier and very different cultural and material adaptation (electronic supplementary material, figure S3). The algorithm also associates additional unit-levels at other sites with CP 14 that, consistent with the hypothesis of these coming from an earlier culture, occur largely at the greatest depths within EUs, as in sites 631 EU 1, 636 EU 4, 642 EU 4. These deeper unit-levels may be associated with a Neolithic-to-Iron Age cultural transition occurring about 2 500 years ago. Exceptions occur within two sites, 744 and 755 (both located near site 756), with unit-levels associated with CP 14 at near-surface unit-levels, indicating that either these EUs may be partially Neolithic or that taphonomic processes introduced deposits containing Neolithic sherds into later stratigraphic contexts.

### (c) Uneven transitions

Despite the high correlation among successive depths, many sites also exhibit repeated switching between two CPs, as exemplified by site 602 EU 2, site 635 EU 3, and site 642 EU 2. These ‘switchbacks’ may indicate slow, inconsistently progressive cultural shifts, potsherd movement between levels, or aggregation of culturally distinct EUs. Simulation results described below indicate that mixture due to potsherd migration would need to be extensive to entirely explain these effects. Consideration of the culture painting for the RU unit-level aggregation suggests that some fraction of these events may be attributed to aggregation effects (see, for instance, site 642 in the electronic supplementary material, figure S8). The remaining events appear to represent shifts within the data, suggesting that the algorithm is in certain cases sufficient to resolve cultural dynamics at previously inaccessible levels of granularity. These will need to be further studied by broader comparisons with other, less common data types, including stone tools and clay figurines, but those comparisons are beyond the boundaries of the present analysis.

### (d) RCDs at site 642 indicate ancient segregation

The culture painting for site 642, the most complex and extensive excavation in the region, also exhibits the most complex patterning, with evidence of significant internal differention among EUs. This site is comparable in size to modern farming villages in the region, and such villages often include ethnically segregated neighbourhoods. The algorithm indicates that the units can be organized into five broad groups: (i) EUs 1, 1B, 1C that exhibit CPs 4 and 9 near the surface, transitioning to CPs 3 and 8 at greater depths, (ii) EU 2 has repeated successive shifts between CP1 and CP2 with progression to CP 10 at 2.2 m, (iii) EUs 3 and 4 associate with CP 11 from the surface down to approximately 1.4 m, then change largely into CP 8, (iv) EU 5 shows CP 8 changing to CP 9 and then CP 1, and (v) EUs 6–14 show CPs 3 and 11, with introgressions from CP 8 and significant posterior uncertainty consistent with their low number of counts. Site 631 shows similar but less eleborate segregation patterns between EUs 1 and 2 and EU 3. These stable but distinctive patterns within each group are consistent with the previously hypothesized proposition that modern ethnic segregation within Lake Chad Basin villages is inherited from deeper historical divisions [11].

Ancient segregation patterns are also supported by integrating RCD information with the culture painting for the subset of EUs with that information, as shown in figure 4. EUs 3 and 6 exhibit distinct CPs from both EUs 1 and 2 and from each other for depths dated to approximately 1 000 years before the present, with the culture painting showing consistent CPs for neighbouring unit-levels. EUs 1 and 2 also have largely distinct CPs patterns from the surface down to the dated depths, consistent with the continuity of segregation patterns into the historical period. However, CP2's substantial but variable presence in both EUs suggests more complex cultural dynamics between these two locations than strict segregation.

### (e) Simulations

To ensure that model inference is reliable, an additional simulation study was performed to understand how the total number of potsherds, the number of decoration types, the number of unit-levels, the degree of autocorrelation among levels, and the amount of potsherd mixing across levels contribute to algorithm performance. Data were simulated as follows. The number of CPs was fixed to five and the number of unit-levels was fixed to 100 and organized into five sites with 20 unit-levels. The remaining parameters were varied for each simulation. Unit-levels were assigned randomly to a CP with a probability *ρ* of switching cultures between successive unit-levels. For each CP, a set of parameters were randomly generated to determine the external decoration distribution, with the concentration parameters sampled independently from an exponential distribution with mean one. Conditional upon the total number of counts, decoration data were then generated for each unit-level. To simulate potsherd mixing, a fraction *f* of each unit was swapped with neighbouring unit-levels. A complete description of the parameters used is given in the electronic supplementary material, table S1. Each set of parameters was repeated for 10 iterations. The simulation study only covers the regime similar to the one encountered in the real data, where the number of decoration types is most often significantly smaller than the total number of counts. Given a fixed number of potsherds, extremely high numbers of decoration types might lead to reduced inferrential power but we did not simulate those conditions for these cases explicitly.

As shown in the electronic supplementary material, figure S9, we find that the algorithm's performance depends largely on the total number of potsherds, the number of decoration types, and, to a lesser degree, the amount of mixing between unit-levels, as measured by three metrics: the Kullback–Liebler divergence between the maximum-likelihood posterior sample and the simulated CP assignment; the mean correlation between *A*_{c} for inferred and simulated CPs; and the inferred number of CPs. Performance is broadly consistent for each metric, with little improvement in performance when there are more than 15 decoration types and 250 potsherds within a unit-level. Substantially below these levels inference can be highly variable. The amount of mixing only notably affects inference at very high rates (50% of potsherds) and with low amounts of unit-level correlations. As would be expected from the absence of correlation modelling, performance is insensitive to the degree of autocorrelation among sites when mixture is not considered. Within the Mandara dataset, we have 13 decoration types, with a median number of counts of 325 (7 and 1 567 for 5% and 95% values, respectively). Conditional upon having low mixture (less than or equal to 10%), simulations indicate 90% or better correlation between simulated and inferred CPs for unit-levels with greater than 250 total counts and 7, 15, or 25 decoration types. These criteria are met by 370 out of 516 unit-levels in the Mandara dataset.

## 4. Discussion

The model presented provides a new approach to modelling distributions of ancient artefacts, admitting quantitative comparison about cultural affiliations, their successional dynamics, and the distribution of their material production. Inferred from potsherd data collected at a set of sites ranging across a culturally complex region, the model provides insight that is both consistent with current research and provides significant extensions to current practice, demarcating cultural shifts as well as the uncertainty in these estimates. Advantages over statistical techniques previously used to analyse ceramic data include the avoidance of *a priori* assumptions about the underlying number of model components and their parameters, as well as the ability of the algorithm to resolve very fine-grained variation in CPs.

As we have noted above, CPs are complex analytical entities that are probably derived from a number of interacting sociocultural, environmental, and taphonomic processes. The cultural switching between CPs that we detect should not be thought of as necessarily involving switching between ‘ethnic groups’ in any simplistic sense; the model is simply indicating that the most parsimonious explanation of the data are to treat the succesive CPs as distinct. The shifting potsherd distribution may be attributable to cultural processes within particular ethnic groups, the transcending or subversion of ethnic boundaries, or changes in available resources. This is certainly the case for modern groups in the region [11].

This analysis indicates that aspects of cultural patterning seen today (especially in terms of the spatial distribution of ethnolinguistic groups, and cultural segregation within nucleated villages) may have deep-time antecedents extending more than 1 000 years back into the Iron Age. This is consistent with overall relationships between modern and Iron Age ceramic assemblages in the region. Furthermore, this provides support for a qualitative assessment of differences between Neolithic and Iron Age ceramic assemblages, and for the existence of transitional assemblages in the mid-third-millennium BP. This approach will likely be useful in exploring other cultural processes in the region, including the Iron Age occupation of the Mandara Mountains themselves. Decoration data can also be compared to other regions of Africa in order to examine cultural dynamics over much larger spatial scales.

This modelling approach is not limited to potsherd decoration and may be usefully extended to a variety of artefactual data. In particular, the model only leverages one form of classification, while the majority of potsherds have additional categorization based on vessel type, potsherd type, and potsherd location within the vessel. These data provide complementary information to external decoration so model extensions that include them may provide increased resolution of cultural shifts. At the other extreme, the insistence on high-frequency decoration types in this model also excludes exotic potsherds and decoration types that researchers have traditionally used to infer long-range trade patterns, such as the sgraffito decoration type often held to indicate Kanuri presence or contact [39]. The Bayesian formulation of the model holds out the possibility that these data and other non-ceramic forms of data—notably linguistic, genetic, mineral trace information, and other artefact types (lithics, metals)—could be integrated to provide rich frameworks for data analysis and hypothesis testing.

The integration with RCDs presents a point of particular interest to archaeology, as it connects physical depths with absolute dates. As shown in figure 3, the model's inference provides a sequential understanding of cultural patterning within each EU that is not always consistent within sites. Extending the model's statistical scope to include these data and map all unit-levels onto a uniform timescale would further enhance the portrait of ancient taphonic and cultural dynamics. This would also provide a new method for understanding the interaction of settlements, soil types, and deposition rates, potentially between and within sites. However, this extension presents signficant statistical challenges since bringing conflicting RCDs and unit-levels that lack RCDs into correspondence has no ready solution.

The analysis undertaken here deliberately does not model correlations among neighbouring unit-levels or across EUs within the same site. This is partially to understand the sensitivity of the model to detect distributional shifts, but also to avoid the perils of overfitting. However, modelling these correlations provides the promise of ‘borrowing strength’ across EUs to impute CPs, where data are scarce, as is commonly possible in Bayesian analysis [40]. The specification of infinite hidden Markov models, particularly allowing for variable dwelling times, may prove a fruitful extension [41]. The generation of parameters for the DM distribution may also be extended by modelling how more recent cultures arise from older ones according to, for instance, a hierarchical Dirichlet process, or by providing additional parameteric flexibility afforded by a logistic multivariate normal distribution [42].

An operational question raised by this study is how researchers should aggregate potsherd data for analysis by the culture painting algorithm. As noted in the Introduction, in the Mandara mountain excavation most sites contain multiple EUs that in turn often contain multiple RUs. There is no intrinsic feature of the model other than the number of counts for each unit-level that suggest at what level the data should be combined to provide the most relevant archaeological understanding. The model provides largely similar results when aggregated at the EU or RU level, but with signficantly more uncertainly at the lower level of aggegation. However, apparent contradictions in the RCD dating at the EU-level—where lesser depths yield earlier dates—are sometimes resolved by disaggregating to the RU level (electronic supplementary material, figure S8). Aggregation at the site level leads to fewer inferred CPs and markedly less discrimination among the sites, showing the pitfalls of overaggregation (electronic supplementary material, figure S6). However, these observations may be specific to the context of the Mandara mountains and general guidance as to how to navigate this tension will only be possible through corroboration with other artefact types and the method's application to a range of archaeological contexts.

## Authors' contributions

J.D.O. designed the study, implemented the statistical methods, and co-wrote the manuscript. K.L. implemented the visualisations and edited the manuscript. S.M. assisted in the study design and co-wrote the manuscript.

## Competing interests

We declare we have no competing interests.

## Funding

We received no funding for this study.

## Acknowledgements

The authors are grateful to Gabrielle Grabin, Ana Lagunez, and Ryan Knauss for careful reviews and thoughtful discussion of the manuscript.

- Received November 23, 2015.
- Accepted February 24, 2016.

- © 2016 The Author(s)

Published by the Royal Society. All rights reserved.