Technical and Methodological Comments on McLaughlin et al. (2018)

McLaughlin et al. (2018) argue for similar trends in the medieval Irish historical record and the archaeological radiocarbon record. Part of their results are due to an ad-hoc bandwidth being used to calculate the kernel density estimates (KDEs). This contribution looks at using a data-driven bandwidth to re-calculate the KDEs and look at the first derivative of the KDEs. The results here indicate the radiocarbon record declines much sooner than the early 9th century and not recovering again until the late eleventh century. Comments are also noted on the Irish annals and the approach evidenced in the archaeological literature, for at least one region, on the use of radiocarbon dates.


Introduction
While McLaughlin et al. (2018;hereafter ML) is an ambitious and interesting approach that works at linking the archaeological data to the historical annals for medieval Ireland, there are some methodological and technical shortcomings that need to be noted. The purpose of this research note is to highlight some of the methodological and technical issues with ML and to provide a re-analysis of the data.

Further Considerations: Radiocarbon Dates
The radiocarbon record depends on what survives that can be dated, and what archaeologists determine need to be dated. ML (14) argue that the size of the radiocarbon database renders null any bias, though they fail to fully acknowledge excavator and research bias. From a review of the Ulster Journal of Archaeology from 2000 through 2016, 24 sites with a High Medieval (after 1169) component were excavated. Only 13 of the 24 sites with a High Medieval component were dated by radiocarbon dating; in many cases the sites or their particular component were dated primarily by pottery or other finds, and a single radiocarbon date is often only provided as supplemental information. For Early Medieval sites over the same time, out of 17 sites, 13 were radiocarbon dated. At many of these sites, 1 While there are those who disagree with McCarthy (see for example Charles-Edwards 2009/2010), even rejecting his reconstruction of the Annals of Ulster, there is an earlier set of scholars which advanced the Annals of Ulster were a compilation of earlier annals (Smyth 1972). 5 multiple radiocarbon dates of the Early Medieval component were taken. 2 An examination of the potential impact of radiocarbon dating more sites with a High Medieval period component will be discussed in a following section.
Also of concern is the quality of the dataset under consideration. How thorough has been the search for dates from Irish sites, or has it depended only on what was in Chapple (2015)? ML only lists ten dates from the site of Ballyhanna, Co. Donegal in their database. Going to the original publication for Ballyhanna (McKenzie and Murphy 2011), one will find that there are 34 radiocarbon dates for the site that date between 1800 and 600 radiocarbon years bp. 3

Kernel Density Estimates
Kernel density estimates (KDEs) have recently been adopted by those working with radiocarbon dates to reduce spurious peaks caused by the calibration curve (Bronk Ramsey 2017;Brown 2015;. Average or composite estimates of multiple KDEs are used to reduce calibration artifacts and potentially are a better reflection of the true growth rate (Brown 2017). Individual KDEs based on a single sample do not conflate process variation and chronological uncertainty, but when the samples are combined via bootstrapping (or any other relevant method), the resultant combined KDE does mix chronological uncertainty and process variation. This means KDEs can be useful when employing an "integral approach" and less so for a "point-wise approach" (Carleton and Groucutt 2021, 630-34). 4 One oversight in the use of KDEs by some archaeologists is the lack of a mathematical or statistical methodology in calculating the bandwidth which is then used to calculate the KDE. 5 This issue received attention in Bronk Ramsey (2017), and a methodology for dealing with the bandwidth issue is provided in the 2 Crema (2021c, 5) notes the bias in using other dating methods when dealing with historical periods. 3 Twenty-four additional samples would increase the burial data set for Ireland by almost 5%. My interest is in the archaeology and history of the province of Ulster. ML's database contains only 470 radiocarbon dates from Ulster between the years 1800 and 600 radiocarbon years bp. For the same radiocarbon time period, I have assembled a database containing 510 radiocarbon dates that were published between 1960 and 2016. 4 The "integral approach" focuses on the temporal distribution of the data set and utilizes areas of summed probability distributions (SPDs) or KDEs across intervals of time. The "point-wise approach" focuses on the number of events dated to a particular time. See Carleton and Groucutt (2021) for a more in-depth discussion. 5 For example, an ad-hoc bandwidth is used in Brown and Ames 2019, Capuzzo et al. 2020, andRiris 2019. OxCal program (Bronk Ramsey 2009). Bandwidth selection is more recently discussed in Crema (2021c).
In the case of ML, an ad-hoc bandwidth was used to calculate the KDEs. Sheather (2004) noted the bandwidth is the most crucial element in determining the shape of the KDE. Ideally, a properly chosen bandwidth is one that minimizes the mean integrated squared error (MISE) (Hazelton 1996, Baxter et al. 1997, Jones et al. 1998, Sheather 2004. When the bandwidth is too small, it will under-smooth the KDE, create spurious peaks, and contain noise (Loader 1999, Sheather 2004. When the bandwidth is too large, over-smoothing will occur and features in the KDE will be obscured (Sheather 2004).
For the results discussed in this contribution, a data driven bandwidth as proposed by Sheather and Jones (1991;cf. Sheather 2004) is calculated from each bootstrap sample of the calibrated radiocarbon dates. While the Sheather-Jones plug-in method is known to over-smooth complex densities (Loader 1999), its overall performance provides an optimal balance between bias and variance (Jones et al. 1996).
The first derivative of a kernel density estimate is the instantaneous rate of change at any given point in time. While ML discusses trends in the data, there is no discussion in the rates of change in the KDE. The rates of change are important when discussing downward trends and stability. In the re-analysis below, first derivative plots will be presented.
The 95% confidence intervals (CI) and the median are calculated from the bootstrap samples and are illustrated in both the KDE and first derivative plots. As a reminder, the 95% CI is an interval estimate that we can be 95% confident, out of all intervals, contains the parameter of interest (Gardner and Altman 1986;Sim and Reid 1999). It also bears noting that KDE confidence intervals are point estimates (Chen 2017:164, 167-69).

Re-analysis of the data set
The re-analysis presented below examines: 1) all the radiocarbon dates from Ireland in ML; 2) just the burial radiocarbon dates from Ireland; 3) radiocarbon dates of cereals from a variety of Irish contexts; and 4) the settlement debris. Bootstrapped KDEs and first derivatives of the KDEs were generated from the calibrated radiocarbon dates.
All statistical operations were done in the open-source package R (R Core Team 2020). Uncalibrated radiocarbon dates were first calibrated using Bchron (Haslett and Parnell 2008), and then calibrated probability distributions were sampled 1000 times. The bandwidth of each sample was calculated using the Sheather-Jones plug-in method as implemented by the hsj command in the sm package (Bowman and Azzalini 2018). KDEs and the first derivative of the KDE were calculated using the kedd (Guidoum 2015) library. Global modes are calculated using the hdr command from the hdrcde library (Hyndman, Einbeck and Wand 2021). The overlapping coefficient (OVL) was calculated with the sfsmisc library (Seminar fuer Statistick 2021). Figure 1 contains plots of the re-analysis of all the radiocarbon dates for Ireland. The histogram of the calculated Sheather-Jones plug-in bandwidths from the 1000 bootstrap samples is in the left panel. The mean bandwidth for all the data is equal to 42.68 years, with a 95% CI ranging from 37.56 to 47.24 years. The right panel in Figure One is a comparison of the two KDEs. The 95% CIs are represented by a dotted line, while the solid line represents the median value of the KDE curve. The black curves are calculated with a 30-year bandwidth and the red curves are calculated with a Sheather-Jones plug-in bandwidth. What is noticeable in both figures is that the confidence interval--which is a measure of variability--is much narrower with the Sheather-Jones plug-in bandwidth.  Simulated-15% n/a n/a 659 659 Simulated-30% n/a n/a 1255 1258 For the KDE calculated with the 30-year bandwidth (KDE30), based on the 95% CI, the KDE peaks sometime between 669-673 CE ( Table 1). The KDE calculated using the Sheather-Jones plug-in (KDESJ) peaks between 685-688 CE. After both KDEs peak, they begin to decline. Looking at the first derivatives of the KDE30, the 95% CI becomes negative between 654-855 CE ( Figure 2). Thereafter, the 95% CI contains negative rates of change. For the KDESJ, the rate of change becomes negative between 666-826 CE. Likewise, the 95% CI for this first derivative also contains negative values after 826 CE.
For the entire dataset for Ireland, the KDEs and first derivatives, despite being calculated with different bandwidth approaches, have similar results. Both KDEs contain a global mode dating to the late 7th century and the densities decrease thereafter. The 95% CI for the first derivative of the KDEs have negative rates of change between the second half of the 7th century and the first half of the 9th century.

Figure 2.
First derivative (with 95% CIs) for the two KDEs. The rate of change becomes negative sometime between the mid-7th to 9th centuries CE. Note the width of the CIs; the one calculated with a Sheather-Jones plug-in bandwidth is much narrower.
Figures 3 and 4 plot the results for the radiocarbon dates from burials. The mean bandwidth is equal to 52.79 years, with a 95% CI ranging from 43.44 to 60.47 years. The KDE30 has a global mode between 625-638 CE ( Table One). The KDESJ peaks between 642-649. The first derivative of the KDE30 becomes negative between 608-656 CE (Figure 4). The 95% CI contains both positive and negative values thereafter. For the KDESJ, the 95% CI for the first derivative becomes negative sometime between 624-662 CE. Between 1049-1163 CE, the rate of change becomes positive again and stays that way until sometime between 1216-1262 CE. Both KDEs point to a downturn in the 7th century.     The re-analysis of the settlement debris is in Figure 7. The mean of the Sheather-Jones plug-in bandwidth is 80.64 years, and the 95% CI is 66.75 to 89.93 years. Once again, the median KDEs illustrate the issue with undersmoothing using a 30-year bandwidth. The KDE30 peaks between 836-841 CE. The first derivative of KDE30, based on the 95% CI, becomes negative between 657 and 1264 CE ( Figure  8). Thereafter the 95% CI covers both positive and negative values. For the KDESJ, it peaks between 683-737 CE. Negative rates initially occur sometime between 763-861 CE in the first derivative of KDESJ. The 95% CI stays negative after that. The first derivative of KDE30 becomes negative sometime between 657-1264.

High Medieval period radiocarbon dates
As noted above, the lack of dating High Medieval sites by radiocarbon dating can influence the results. For the data set under consideration, 2331 out of 2636 radiocarbon dates pre-date the start of the High Medieval period (circa 880 radiocarbon years bp or 1169 CE).
To look at the effect of dating more samples from the High Medieval period, two simulated data sets were created from the original radiocarbon data. The first simulated data set is created by randomly sampling with replacement the raw radiocarbon dates and standard deviations to create a new set of radiocarbon determinations such that 15% of the radiocarbon dates post-date the start of the High Medieval period. The second simulated data set was created in similar fashion so that 791 out of 2636 radiocarbon dates, or 30%, post-date the start of the High Medieval period. Figure 9 contains the plots of the KDEs created by using the original and simulated radiocarbon data sets. The mean bandwidth for the 15% data set is 37.92, with a 95% CI ranging from 33.72 to 42.07. The mean bandwidth for the 30% data set is 26.63, with a 95% CI ranging from 24.85 to 28.40. Figure 9. KDEs for the simulated data sets versus the original data set. All KDEs were calculated using a Sheather-Jones plug-in bandwidth. As illustrated, the increasing number of High Medieval period dates creates a new mode in the KDE. A simulated data set with 35% High Medieval period dates has a significantly different global mode.
The differences are readily apparent between the three KDEs. As the radiocarbon dates increase for the High Medieval period, a new mode appears in the 12 th to 13 th centuries CE. The global mode for the 15% data set is CE 659 (see Table One). Despite having an OVL ranging from 0.958 to 0.962, there are statistically significant differences between the KDEs of the original and 15% data set since their confidence intervals do not overlap between circa 1150-1350 CE. 6 The global mode for the 30% data set is 1255-1259 CE. Comparing the 30%data set to the original, the OVL ranges from 0.809 to 0.839, and the two 95% CIs do not overlap between circa 500-950 CE and circa 1150-1450 CE.
These two simulations for the High Medieval period illustrate the effect of biases in sampling and dating. Price et al. (2021, 2) note that even as the number of observations increase for a KDE or a summed probability distribution (SPD), they will not converge on the true distribution. Sampling bias can exacerbate this issue even further.
One is justified in asking if there is suitable material to date from the High Medieval period in Ireland. The answer is "yes". Some examples are, but not limited to: a) the faunal remains from the multi-period site of Aghavea, Co.

Summary
The main points concerning ML's use of data are: 1. While the Annals of Ulster have time depth, there are other suitable annals that could be used in a study like that which ML offer. The Annals of Ulster were subject to political tampering as the other annals. 2. While ML note some shortcomings with the radiocarbon record, they fail to note that at least in Ulster there is a propensity to date Early Medieval sites with multiple radiocarbon dates, while High Medieval period sites are less commonly radiocarbon dated. The lack of dates in this later phase, as demonstrated with two simulations, can impact the modes of the KDE. 3. The 30-year bandwidth in all cases under-smooths the dataset and creates spurious peaks in the KDE. This is particularly evident in the cereal and settlement debris KDEs. The Sheather-Jones plug-in method minimizes the variance of the KDE, and thus, the 95% CI for the resulting KDE and its first derivative are less than those that were created with an ad-hoc, 30-year bandwidth. 4. Except for the settlement debris radiocarbon dates, the radiocarbon datasets re-analyzed using a Sheather-Jones plug-in bandwidth have negative rates starting in the mid-7th to early 8 th centuries CE. A negative rate implies the curve is now concave and heading towards a minimum, and there is a negative growth rate. ML's argument for stability ML (13) appears to be based solely on the shape of the KDE, overlooking the rate of change as evidenced in the first derivative of the KDEs. This decline is also earlier than stated in ML (11). 5. While ML (13) do not define stability or decrease in their Table 1, I interpret the period between 700 and 825 as unstable since the 95% CI of the first derivative of the KDESJ goes from being positive to negative. This would indicate that there is a declining growth rate. 6. Depending on the radiocarbon dataset examined, the first derivate of KDESJ becomes positive sometime in the late 11 th century to early 12 th century and continues until the mid-13 th century. From a radiocarbon perspective this indicates an upswing in socio-cultural activity. This could possibly be due to the Anglo-Norman/Anglo-Welsh colonization of Ireland.
Prehistorians and others working with differences and uncertainty at the centennial scale may question if the differences at the decadal level or below are truly important. This author would say "yes". Bootstrapped KDEs account for sampling error, chronological uncertainty, and calibration effects (Carleton and Groucutt 2021, 136;Crema 2021c, 10). By utilizing the historical record and palaeoclimate/palaeoenvironment proxies within an "integral approach", one can formulate hypotheses on climate, historical events, and their potential effect on socio-economic development in Ireland during the Early and High Medieval periods.
For example, looking at the global modes and their subsequent decreases for all the radiocarbon dates from Ireland, the decade between 669-678 CE, which encompasses the global mode of KDE30, reveals that Ireland had reconstructed Palmer Drought Severity Index (PDSI) values indicative of a moist climate (Cook et al. 2015). The Annals of Ulster (Airt and Mac Niocaill, 1983), in addition to containing records of dynastic conflicts, note a famine (670 CE), snowfall (670 CE) and later an abundance of mast in this decade (672 CE). The decade between 685 and 694, which contains the global mode for KDESJ, is marked by: a) a decade long mild drought in Ireland (Cook et al. 2015); b) dynastic conflict 7 ; c) an earthquake and windstorm is noted as occurring in 685 CE in the Annals of Ulster, Annals of Tigernach (Niocaill 2011) and the Chronicon Scotorum (Niocaill 2010).
It also needs to be noted that this data set could be re-analyzed by other approaches and methods which may possibly provide alternative results and interpretations. Crema (2021c) divides the possible approaches down into three different groups: 1. reconstructive approaches 2. null-hypothesis significance testing (NHST) 3. model-fitting approaches ML and this study are characteristic of the reconstructive approach. NHST approaches, as deployed in rcarbon (Crema and Bevan 2020), allow one to test specific hypotheses through the use of Monte Carlo simulation and permutation tests. Model-fitting approaches have been more recently developed, and some examples of programs implementing this approach are ADMUR (Timpson 2020;Timpson et al. 2021), bayesdem (Price et al. 2021), and nimbleCarbon (Crema 2021a(Crema , 2021b. Model-fitting approaches, however, can be computationally intensive, particularly for large data sets.
Which of these methods is best depends on the question being asked. The reconstructive approach, as used here, can provide a quick preliminary, qualitative assessment of differences along the time-axis. While explanations are proffered in ML, no formal hypothesis testing nor modelling occurred. The reconstructive approach only provides a starting point for other analyses.
NHST approaches allow specific hypotheses to be tested-such as whether the distribution of dates conform or deviate from a specific population growth model or whether two sets of dates have similarly shaped KDEs or SPDs. For example, the burial and settlement data presented in ML could be tested to see if the two data sets truly come from similar distributions and have a similar pattern of growth and decline. The NHST approach is more robust but has limitations regarding p-values (Crema 2021c). Like any null hypothesis test, large sample sizes can provide a low p-value-suggesting significance-but the p-value does not assess the effect size nor the probability of the model being correct. 8 Model-fitting approaches allow one to deal with the chronological uncertainty and/or the sampling errors associated with radiocarbon dates. These approaches allow for multiple models to be tested and determine which model fits best. For the Irish radiocarbon data set, for example, nimbleCarbon could be used to determine the changepoint in the growth rate of the SPD of the calibrated radiocarbon dates. The calendar year(s) of the changepoint could then be compared to the dates for when the first derivative of the KDEs become negative.