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