Assessing the effects of age-at-death assignment and allocation on the interpretation of past sheep/goat management practices

This is the Material and Method section taken from the main manuscript (Gillis et al.). 1 Datasets This study uses two archaeological datasets to explore methodological variations within the Payne age-at-death system. Firstly, an unusually large sample of raw caprine mandibular wear data from the Turkish prehistoric site of Çatalhöyük was used to test the impact of inclusion criteria and wear-stage assignment methods on age at death results at a single site. This dataset combines published results from the Neolithic East Mound (Russell et al. 2013; Russell et al. 2014; Twiss et al. 2021) with new data from the Chalcolithic West Mound (Orton in press). Secondly, sheep and goat teeth and mandibular data assembled from 13 early Neolithic sites in Central Europe and the northwestern Mediterranean dating to the 6th-5th millennium BC, were used to compare different allocation methods and visualisation options when comparing results from more fragmented assemblages. This dataset overlaps with those used in previous papers on caprine age-at-death (Gerbault et al. 2016; Timpson et al. 2018). See Table 1: Sheep and goat mandibular data assembled from 13 early Neolithic sites from Central Europe and the northwestern Mediterranean, plus Neolithic and Chalcolithic Çatalhöyük in central Anatolia. ALP: Alförd linear pottery; LBK: Linear Bandkeramik; PN: Pottery Neolithic; EC1: Early Chalcolithic I. 2 Age-at-death determination methodology Age-at-death was determined from teeth (dP4/P4, M1-3) using eruption, replacement and wear stage (Payne 1973). The analysis was carried out by DO (Çatalhöyük) and REG (Neolithic Euroepan sites). Isolated M1 and M2 from heavily fragmented Neolithic assemblages were distinguished both by visual comparison and metrically using the height of the crown at the enamel junction and breadth of the tooth at the base of the crown (Vigne unpublished) similar approach as Beasley et al. (1993) with cattle. For the Çatalhöyük dataset, four approaches to assigning Payne stages were applied, from the most conservative to the most liberal: (i). Directly staged only: specimens that could be placed directly into a specific age category based on Payne’s core criteria (equivalent to ‘raw count’ in Payne 1973, Table 1). (ii). Core criteria (“minimal”): all specimens that could be assigned to one, two, or three stages based on Payne’s core criteria. (iii). Assuming dP4/M3 exclusivity (“conservative” (‘cons’)): all specimens that could be assigned to one, two, or three stages based on Payne’s core criteria with the added assumption that dP4 and M3 are never in wear at the same time. (iv). Using Payne’s tables (“extrapolated”): all specimens that could be assigned to one, two, or three stages either directly or by extrapolation from preserved teeth using Payne’s (1973, Figures 11-13). This analysis was repeated firstly using all mandibles and secondly only using those with a dP4 or M3 present (following Payne’s criterion), for a total of eight variants. All eight of these were used for an initial comparison of results, and the results of three (approaches 2, 3, and 4 using the dP4/M3 inclusion criterion) were carried forward for downstream analysis (see below). For the Neolithic multi-site dataset, age stages were assigned using the full method from Payne (1973), including extrapolation from preserved teeth (equivalent to assignment approach 4 above), but without applying the dP4/M3 inclusion criterion, as no mandibles had in wear dP4 at the same time as M3. All datasets were initially summarised as minimum number of individual counts (MNI) without allocation to specific classes (c.f. Timpson et al. 2018), capturing that some teeth/mandibles can only be assigned to several age classes (e.g. ABC; see Supplementary Table 1). In all cases, only the best-represented side of the body was included to avoid double-counting. 3 Allocation methods The three allocation methods described above were then applied to the Neolithic multi-site dataset to tackle the problems of isolated teeth and fragmented mandibles that could not be assigned to specific wear stages. (i). Payne: allocating specimens across possible age classes in proportion to the relative frequency of definitively-staged mandibles in each class, i.e. “proportionally allocated” as set out by Payne (1973); (ii). Flat: allocating specimens equally across all possible age classes (A=0.5, B=0.5); For example, if a tooth/fragmented mandible was summarised as belonging to stage “ABC”, then in the flat allocation process, age classes A, B and C would receive an equal proportion of that element, i.e. (⅓=) 0.333 each. This was referred to as “direct allocation” by Arbuckle 2006. (iii). Duration-based: allocating specimens across possible age classes in proportion to the length (in months) of those classes (e.g. A=0.2, B=0.8). For the same example above, the duration-based allocation process would instead consider that since age class C is the longest (in months) of the three, this age class would receive more of that mandible than the other two younger and shorter age classes. More specifically, following the example above, considering the total duration of stage “ABC” is 12 months, where A, B and C last two, four and six months, respectively, the duration method would allocate (2/12/1=) 0.166, (4/12/1=) 0.333 and (6/12/1=) 0.5 of an isolated tooth/fragmented mandibles to age classes A, B and C, respectively. The same three allocation methods were then applied to three variants of the Çatalhöyük dataset, namely assignment approaches 2, 3, and 4 above, applying the dP4/M3 inclusion criterion, in order to explore interactions between stage assignment approaches and allocation methods. 4 Implementation We believe the strength of the approach we present lies in the scripting of the three allocation methods (and other downstream analyses) in a single programming language script (R). Scripting these methods makes allocation of teeth to age classes transparent, allows discussions on the decisions involved in the allocation process to be more explicit, and will consequently promote improvements and standardisation. Both the stage assignment process for Çatalhöyük and the allocation processes for the multi-site dataset were implemented in R v. 4.3.2 (R Core Team 2023) using the data.table package (Barrett et al. 2025) and custom functions developed by the authors. R scripts and data are available in the open access repository. 5 Age-at-death profile construction 5.1. Monte-Carlo approach In order to visualise differences in age distributions due to assignment approaches applied to the Çatalhöyük dataset (i.e. before any additional divergence introduced by allocation process), a Monte Carlo approach was used to construct survivorship curves for each assignment approach detailed above (i.e. directly staged, minimal, conservative, extrapolated). Each specimen was assigned a random age in days within the limits suggested by its wear stage; for example, a specimen staged at “BC” – corresponding to a suggested age range of 2-12 months – might be assigned anything from 62 to 365 days. This process was repeated 1000 times and the results from each repeat plotted as survivorship curves at low opacity on the same axes, collectively visualising the structure of the age data and the associated uncertainty. 5.2. Mortality curves and profiles, and Survivorship curves Here, age-at-death date is represented in three ways: best fit Gamma distributions as continuous function of age (Timpson et al. 2018), histograms (Gerbault et al. 2016) and survivorship curves (Price et al. 2016). Timpson et al. (2018) proposed to compute the probability of the observed data under every possible arrangement of teeth/mandibles observed at a site, by applying a multinomial distribution (AB=1). Teeth counts in Payne’s age classes A to I with multi-class age assignments were used (Supplementary Table 1). Then, maximum likelihood estimates (MLE) of the Gamma distribution mean and shape parameters were obtained with the function gammaMLparameters from the R package GammaModel (Timpson et al. 2018). These joint Gamma shape and mean parameters were finally used to represent the corresponding mortality curve as the best fit Gamma distribution as a continuous function of age. This was performed for each archaeological site of the 13 Neolithic multi-site dataset, (Figure 3A and Supplementary Figure A1), but not for the Çatalhöyük dataset, due to computing-time and -memory required by this method for larger sample sizes (Table 1). We also represent mortality profiles as histograms, one of the traditional means of visually depicting age-at-death data. Histograms have been constructed following the approach presented in Gerbault et al. (2016), except here nine separate age classes were considered instead of seven. Briefly, this Bayesian approach incorporates a Dirichlet distribution to simulate multiple deviates of an observed age-at-death profile. In the context of this study 1,000 deviates were generated per observed age-at-death profile. These deviates were consequently used to create credible intervals for each age class for each observed profile (Figure 3B-3D, and Supplementary Figure 1B-1D). Here, histograms of the age-at-death data and corresponding Dirichlet-derived credible intervals for each age class have been obtained for each of the three allocation methods tested (i.e. Payne, flat, duration). Note, a credible interval in Bayesian statistics is the probability of a parameter to fall within a particular probability. These are synonymous with confidence intervals in frequency statistics. Confidence intervals for each age class have been used in the past, however these are erroneous as they cannot be calculated for each class as it is impossible to know the mean it is difficult due to 1) sample size, often <30; 2) lack of a mean for a specific age class. Mortality profiles as histograms were generated for each of the 13 early Neolithic sites dataset, and for the Çatalhöyük dataset (Figures 3B-3D and 4A-4C, and Supplementary Figures 1B-1D and 2A-2C). Survivorship curves are the inverse cumulative distributions of the age-at-death data and express the relative frequency of specimens that survived a given age class. We used the equation published in Price et al. (2016: equations 1, 2 and 3, and Table 1). The 95% confidence interval represented on the survivorship curves has been obtained from the Dirichlet deviates computed for each site (considering nine age classes) and each allocation method tested. Survivorship curves were represented for each of the 13 early Neolithic sites datasets, and for the Çatalhöyük dataset (Figures 3E and 4D, and Supplementary Figures 1E and 2D). 6 Comparative methods For comparing a set of age-at-death profiles, we applied correspondence analysis on the Dirichlet deviates and fitted Gamma distributions using MCMC (Gerbault et al. 2016; Timpson et al. 2018). Both approaches account for sampling uncertainty, albeit differently. One of the aspects of interest was to visualise how different these two-dimensional representations of sets of age-at-death profiles appeared, but also how the allocation method (Payne, flat or duration) would affect these two-dimensional representations. Both approaches have been described in detail elsewhere, and the procedure undertaken in the context of the current analysis is only briefly outlined below. For the set of 13 early Neolithic sites, a correspondence analysis was performed on the 13,000 Dirichlet deviates and the 13 observed profiles themselves (13,013 rows and nine columns for the nine age classes). This was repeated for each of the three allocation methods (Figure 5A-5C). On each correspondence analysis representation, ellipses were drawn with the stat_ellipse function (ggplot2), setting the options type=”t” and level=0.95 to show the 95% confidence levels for multivariate t-distributions. As an alternative two-dimensional representation of this dataset, the maximum likelihood estimates and confidence intervals of Gamma shape and mean parameters of the 13 profiles using the Markov Chain Monte Carlo algorithm are shown (Timpson et al. 2018). Note this was performed on the 13 profiles with multi-stage assignments (Figure 5D). Ellipses on this two-dimensional representation were drawn with stat_ellipse as mentioned above. For the Çatalhöyük data, the maximum likelihood estimates and confidence intervals of Gamma shape and mean parameters could not be obtained within the timeframe of this analysis, due to the large sample sizes and the amount of computing memory requested by this method, a possible limitation (Timpson et al. 2018). However, after applying the three allocation methods (Payne, flat and duration) discretising the data into nine age classes, correspondence analyses were computed on this data and the corresponding Dirichlet deviates. Ellipses were drawn on the correspondence analysis using the Dirichlet deviates as mentioned above (Figure 6A-6C). 3.7 Statistical assessment of differences in age class assignment between allocation methods In general, in any field of science, interpretation relying upon visual inspection of data patterns often leads to subjective (and sometimes contradictory) story-telling, partly due to the wide range of processes that can affect the data prior to collection and analysis (Gerbault et al. 2014). While it is accepted that tooth allocation to single discrete age classes can artificially skew the pattern drawn (Price et al. 2016, Timpson et al. 2018), this discretisation into age classes facilitates inferring population demographics, without assuming that the exact calendar ages at which skeletal features develop have been conserved across prehistorical and present times (Deniz & Payne 1982; Twiss 2008). Furthermore, various approaches have been devised to make sense of age-at-death data in a given archaeological context, but here our aim was not to decide which one is best. We rather sought to integrate some of those approaches and investigate whether some introduced consistent biases relatively to others. To this end, we explored a statistical means of assessing whether some age classes were over- or under-represented when using any of the three allocation methods (Payne, flat and duration). It could also be possible that no significant differences between age classes can be detected, in which case the use of one or the other should not make much difference to the overall pattern and its description or downstream inferences. We decided to apply linear regression modelling as a statistical test to assess the association between the difference in the number of specimens in age-class between two allocation methods, i.e. duration-Payne, duration-flat and Payne-flat, so that each of the nine age classes were described by three pairwise differences. For example, the difference between duration and flat (duration-flat) for age class B could either be negative, implying more specimens in flat than duration allocation, positive, implying more specimens in duration than flat allocation, or 0, implying no difference between these allocation methods. This example is reported in Duration_Flat or reported as variable Duration_Flat:ageClassB (Supplementary Table 3 and Figure 7). This calculation was done per assemblage and age class for each pairwise comparison of allocation methods (Duration_Payne, Duration_Flat, Payne_Flat). In this linear modelling approach, the outcome variable was the value of the difference, the explanatory variables were the age classes, pairwise comparison of allocation methods and the number of specimens observed (MNI). Details of the model selection procedure and linear regression diagnostics statistics are presented in the Supplementary Tables 2-4 and Supplementary Figure 3). A similar approach was applied to the Çatalhöyük dataset to investigate whether the age determination approaches would have a significant impact on the profiles obtained from the three different allocation methods (Supplementary Figures 4 and 5, Supplementary Tables 5 and 6).