A Bayesian multivariate approach to comparing lead isotope ratios

Abstract

Lead isotope analysis is widely used in archaeological provenance studies, yet comparisons between assemblages remain largely qualitative, relying on visual inspection of bivariate scatter plots. Although Bayesian regression methods have advanced considerably in archaeological applications to radiocarbon chronology and compositional analysis, their extension to lead isotope data has not been systematically explored. This paper addresses this gap by introducing a Bayesian multivariate Student-t regression framework for formally comparing lead isotope ratio distributions (²⁰⁶Pb/²⁰⁴Pb, ²⁰⁷Pb/²⁰⁴Pb, ²⁰⁸Pb/²⁰⁴Pb) between archaeological sites. The model simultaneously accounts for multivariate correlations between isotope ratios, heavy-tailed distributions, and known analytical measurement error. We apply this method to 71 bronzes from Panlongcheng and Zhengzhou Shang City, two key Upper Erligang period sites whose metal provenance relationship has long been debated. The Student-t model is strongly preferred over a Gaussian alternative, indicating strong departures from normality in the data. All three isotope ratios show statistically significant mean differences between sites, with Zhengzhou exhibiting higher mean values. Variances are similar between the two sites. These results provide quantitative support for the hypothesis that Panlongcheng accessed partially different metal sources than Zhengzhou, consistent with the presence of recycled older bronzes at Panlongcheng identified by previous studies. The Bayesian multivariate approach offers a reproducible, probabilistic alternative to visual scatter plot comparison and is readily transferable to other assemblage comparisons in archaeological lead isotope research.

Graphical abstract

Keywords

Lead isotope analysis / Bayesian regression / bronze provenance / Panlongcheng / Zhengzhou Shang City / Erligang period / multivariate Student-t model

1 Introduction

The question of whether regional centres in the early Bronze Age China produced their own ritual bronzes or received them from the political capital has been central to understanding the political economy of the Erligang state (c. 1500–1300 BCE) (Liu, L. and Chen 2003:133). Panlongcheng, located in the Middle Yangtze River valley near modern Wuhan, Hubei Province, was a major Shang settlement and likely the highest-ranking southern outpost during the Erligang period (Zhang 2014). Its relationship with Zhengzhou Shang City, the political capital of the Erligang state in Henan Province, remains a subject of active scholarly debate (Figure 1).

Two principal positions have emerged. One view holds that Panlongcheng served primarily as a node in a metal supply network, receiving finished bronze ritual vessels from Zhengzhou or operating under direct technological supervision from the capital. The other argues that Panlongcheng developed at least partial independence in bronze production, with its own access to metal raw materials and casting capabilities. Resolving this debate requires rigorous comparison of the metal provenance signatures of bronzes from the two sites (Liu, R. et al. 2019; Zhang 2024; Li 2020:72; Chen et al. 2024; Liu, S. et al. 2023).

Multiple lines of evidence have been considered on this question. The analysis of clay cores left in Panlongcheng bronzes shows that they are geochemically different from local deposits and more similar to loess deposits in the north (Nan et al. 2008). This finding supports the argument for importing bronze from the Central Plains. Trace-element-based copper groups have revealed that approximately 30% of Panlongcheng bronzes contain nickel-bearing copper (Groups CG5, CG11) that is rare at Zhengzhou (<8%), suggesting access to different copper sources (Liu, R. et al. 2019). A Kolmogorov–Smirnov test confirmed statistically significant differences in tin content between the two sites, with Panlongcheng vessels containing more tin, likely reflecting proximity to southern tin sources (Liu, R. et al. 2017). The discovery of casting remains at Xiaozui, including crucibles, slag, dross and metal droplets, provides direct physical evidence of local bronze casting activity (Liu, S. et al. 2020; School of History of Wuhan University et al. 2019). Crucible technology also differs between the two sites. Panlongcheng crucibles at the Xiaozui workshop are shallow-bellied vessels made from local iron-rich clay, stylistically distinct from Zhengzhou casting equipment (Liu, S. et al. 2020). However, these scattered lines of evidence could still not rule out that large ritual vessels found in the tombs of Panlongcheng were mostly made at and imported from Zhengzhou.

Lead isotope analysis (LIA) has been a primary tool in this debate. Both sites yield bronzes with highly radiogenic lead (²⁰⁶Pb/²⁰⁴Pb > 19), reflecting a shared source of this distinctive material that characterises many Shang bronzes (Liu, R. and Pollard 2022). However, closer inspection reveals differences in the non-radiogenic lead components. Liu, R. et al. (2019) distinguished two types of common lead: ‘type A’ (²⁰⁶Pb/²⁰⁴Pb ≈ 17.5), more frequent at Panlongcheng, and ‘type B’ (²⁰⁶Pb/²⁰⁴Pb ≈ 16.5), more frequent at Zhengzhou. Critically, Liu, S. et al. (2023) identified a group of less radiogenic lead at the Panlongcheng Xiaozui workshop (²⁰⁶Pb/²⁰⁴Pb ≈ 16.5), arguing that these values reflect recycled older Lower Erligang bronzes originally brought south from the Central Plains. This observation implies that Panlongcheng’s isotopic signature was shaped not only by contemporaneous metal procurement but also by the recycling of earlier material.

Despite these observations, comparisons of lead isotope data between Panlongcheng and Zhengzhou have relied primarily on visual inspection of bivariate scatter plots (e.g., Figure 2). A comprehensive multivariate probabilistic framework for formal assemblage comparison remains lacking. The prevailing scatter-plot approach, while intuitive, is limited in several important respects. First, it is subjective as different researchers may reach different conclusions from the same plots. Second, it cannot quantify the magnitude or uncertainty of differences between assemblages. Third, it treats the three isotope ratios as independent pairs rather than modelling their joint distribution. Fourth, it ignores both the heavy-tailed nature of lead isotope data and analytical measurement error.

These limitations are widely recognised in the broader archaeometric literature. Baxter (1999) demonstrated that multivariate normality is ‘the exception rather than the rule’ for lead isotope data from archaeological contexts. De Ceuster et al. (2023) argued that conventional biplots are infeasible for large databases and called for statistically transparent alternatives. Crema (2025) urged the adoption of measurement-error models beyond chronological applications, noting that ignoring measurement uncertainty can bias parameter estimates. Liu, S. et al. (2025), in a comprehensive review of Chinese bronze provenance studies, explicitly advocated for model-based quantitative approaches to replace visual biplot inspection.

Recent years have seen growing interest in applying Bayesian regression methods to compositional and isotopic data in archaeology. Bayesian approaches offer principled uncertainty quantification, seamless incorporation of measurement error, and natural accommodation of non-normal distributions through flexible likelihood specifications. In ceramic geochemistry, hierarchical Bayesian models have been used to estimate group membership probabilities while accounting for compositional closure and analytical uncertainty (Papageorgiou 2020). In strontium and oxygen isotope studies, Bayesian mixture models have been employed to classify local versus non-local individuals (Scaffidi and Knudson 2020). In chronological modelling, the integration of radiocarbon dates with stratigraphic constraints within a Bayesian framework is now standard practice (Bronk Ramsey 2009; Banks et al. 2019), with recent extensions to questions of technological origins (Rotunno and Crema 2025). In zooarchaeology, Bayesian multilevel models have been applied to mortality profiles (Gerbault et al. 2016), faunal taxonomic abundances (Ragno 2024), and skeletal proportions from commingled remains (Cao et al. 2024). Most recently, Vieri et al. (2025) introduced Bayesian beta regression for compositional variability in metallurgical studies, demonstrating the potential of distributional regression to distinguish technological traditions. However, the application of Bayesian regression to lead isotope data, where the joint multivariate structure, heavy-tailed distributions, and correlated measurement errors present distinctive challenges, has received comparatively little attention. The present study addresses this gap by extending the Bayesian regression paradigm to the multivariate comparison of lead isotope assemblages, introducing a framework specifically designed for the statistical properties of archaeological lead isotope data.

In this paper, we introduce a Bayesian multivariate Student-t regression framework for comparing lead isotope ratio distributions between archaeological assemblages. The model treats all three isotope ratios simultaneously, estimates residual correlations between them, accommodates heavy-tailed distributions through a Student-t likelihood with estimated degrees of freedom, and incorporates known analytical measurement error. We apply this framework to compare 71 bronzes from Panlongcheng (n = 38) and Zhengzhou Shang City (n = 33), both dating to the Upper Erligang period. The approach yields posterior probability distributions for site mean differences, variance ratios, and residual correlations, providing quantitative answers to questions that have previously been addressed only by visual assessment.

2 Archaeological Background

Zhengzhou Shang City, located in modern Zhengzhou, Henan Province, was the political capital of the Erligang state during the early Shang period (c. 1500–1300 BCE; Figure 1) (Liu, L. and Chen 2012: 280-281). The site contains extensive evidence of large-scale bronze production, including the Nanguanwai foundry complex, where crucible slag analysis has revealed three distinct metal supply lines, including copper with common lead from Middle Yangtze sources, highly radiogenic lead introduced primarily with tin and lead, and recycled bronze with low radiogenic signatures (Sun, Z. et al. 2023a).

Panlongcheng, situated approximately 500 km to the south near modern Wuhan in Hubei Province, occupied a strategic position in the Middle Yangtze River valley (Figure 1) (Liu, L. and Chen 2012: 285). During the Erligang period, it was the highest-ranking Shang outpost in the region, with a walled enclosure, elite burials containing bronze ritual vessels, and, as excavations at the Xiaozui locus have confirmed, a bronze casting workshop (Liu, S. et al. 2020; School of History of Wuhan University et al. 2019; Li 2020: 80). The workshop yielded crucibles, slag, copper dross, and metal droplets, providing unambiguous evidence that bronze casting took place at the site.

The metal circulation network connecting these sites and their hinterlands was complex and changed over time. Copper likely originated from mining districts in the Middle Yangtze region, specifically the Jiurui metallogenic district in Jiangxi (Sun, Z. et al. 2023a). Tin sources are less well constrained but multiple sources in both north and south China are potential candidates (Liu, S. et al. 2025). The origin of highly radiogenic lead, a distinctive geochemical signature of many Shang bronzes, remains debated, though recent evidence suggests it was introduced as an additive metal (with tin and lead) rather than as a component of copper ores (Liu, S. et al. 2025; Liu, R. and Pollard 2022). Chen et al. (2024) have synthesised the broader regional picture, describing a network in which metal flow was initially predominantly southward from the Central Plains during the Erligang period but became increasingly bidirectional as southern metallurgy developed greater autonomy in subsequent periods.

The lead isotope data examined in this study are all attributed to the Upper Erligang period. Each sample provides three isotope ratios, ²⁰⁶Pb/²⁰⁴Pb, ²⁰⁷Pb/²⁰⁴Pb, and ²⁰⁸Pb/²⁰⁴Pb. Previous studies have compared these ratios visually using scatter plots, leading to the general impression that the two sites share broadly similar isotopic ranges dominated by highly radiogenic lead. However, Liu, S. et al. (2023) identified a distinct low-ratio component at Xiaozui (²⁰⁶Pb/²⁰⁴Pb ≈ 16.5) that they attributed to recycled Lower Erligang bronzes. This observation, combined with the trace element and alloying evidence summarised above, motivates a formal quantitative comparison of the isotopic distributions of the two assemblages.

3 Methods

3.1 Data

The dataset comprises lead isotope ratios for 71 bronze artefacts, 38 from Panlongcheng and 33 from Zhengzhou Shang City, all dating to the Upper Erligang period. Lead isotope data for Panlongcheng bronzes are taken from the site monograph (Hubei Institute of Archaeology 2001; see also Peng et al. 2001; Sun, S. et al. 2001) and from Liu, S. et al. (2023). Data for Zhengzhou Shang City bronzes are from Tian (2013), Sun, Z. et al. (2023a), and Jin (2008). The full dataset is provided in Table S1. Three response variables are modelled simultaneously. Summary statistics are reported in Table 1. Analytical measurement error is specified as 0.1% relative standard deviation of each individual measurement, a conservative estimate consistent with modern multi-collector inductively coupled plasma mass spectrometry (MC-ICP-MS) and thermal ionisation mass spectrometry (TIMS) precision.

3.2 Statistical Model

We employ a Bayesian multivariate Student-t regression with site as a fixed effect. The model is specified using the brms package (Bürkner 2017), which provides an interface to the Stan probabilistic programming language (Stan Development Team 2023). The model is fitted with family = student() in brms. The multivariate formulation models all three isotope ratios simultaneously, allowing estimation of residual correlations between them.

For each isotope ratio j ∈ {1, 2, 3} corresponding to ²⁰⁶Pb/²⁰⁴Pb, ²⁰⁷Pb/²⁰⁴Pb, and ²⁰⁸Pb/²⁰⁴Pb, and each observation i with site membership s(i), the model is specified as follows.

(1)y_{ij}|{\text{se}}_{ij}\sim \text{Student-}t(\nu ,{\mu }_{s\left[i\right],j},{\tau }_{s\left[i\right],j})

(2){\tau }_{s,j}=\sqrt{{\sigma }_{s,j}^2{\text{+se}}_{ij}^2}

where ν is the degrees-of-freedom parameter shared across all responses and sites, controlling tail heaviness; μ_s,j = β_s,j is the site-specific mean for ratio j; σ_s,j is the site-specific residual standard deviation for ratio j, modelled on the log scale as log(σ_s,j) = γ_s,j; and se_ij is the known analytical standard error for observation i on ratio j. The scale parameter τ_s,j combines the estimated residual variance with the known measurement error in quadrature, so that the total observational variance reflects both analytical uncertainty and genuine between-sample variation.

The three response variables are modelled jointly with a residual correlation structure:

(3)\text{Corr(}{\epsilon }_{i,j},{\epsilon }_{i,k})={\mathbf{R}}_{jk}

where R is a 3 × 3 positive-definite correlation matrix estimated from the data. This multivariate formulation captures the well-known collinearity among lead isotope ratios (Albarède et al. 2012) and permits inferences about between-site differences that respect the joint distribution of all three ratios.

The model uses site-specific sub-models for both the mean (μ, via β) and the log-scale dispersion (log σ, via γ), estimating separate location and scale parameters for each site on each ratio. This parameterisation allows formal comparison of not only central tendency but also variability between assemblages, a distinction that is important for provenance inference but rarely quantified in existing lead isotope studies.

3.3 Priors

Weakly informative priors are specified for all parameters (Table S2):

(4a)\begin{array}{c}{\beta }_{s,j}\sim \text{Normal(0,10)}\end{array}

(4b)\begin{array}{c}{\gamma }_{s,j}\sim \text{Student-}t(3,0,2.5)\end{array}

(4c)\begin{array}{c}\nu \sim \text{Gamma(2,0.1)}\end{array}

(4d)\begin{array}{c}\mathbf{R}\sim \text{LKJCorr}\left(2\right)\end{array}

The Normal (0, 10) priors on site means are effectively non-informative on the scale of lead isotope ratios (typically ranging from 15 to 42). The Student-t (3, 0, 2.5) priors on the log-scale dispersion parameters are weakly informative, regularising extreme variance estimates while permitting a wide range of plausible values (Bürkner 2017). The Gamma (2, 0.1) prior on ν places most prior mass between 2 and 50, accommodating distributions ranging from heavy-tailed to approximately normal. The LKJ (2) prior on the residual correlation matrix mildly favours smaller correlations while remaining broadly permissive (Lewandowski et al. 2009).

3.4 Computation

Posterior samples are obtained via Hamiltonian Monte Carlo (HMC) using Stan, with 4 chains of 4,000 iterations each (2,000 warmup). The target acceptance rate is set to 0.999 and the maximum tree depth to 15 to ensure reliable sampling for this complex model. A fixed random seed (42) ensures reproducibility.

3.5 Model Diagnostics

Convergence is assessed via the potential scale reduction factor (\hat{\mathrm{R}} < 1.01 for all parameters) and effective sample size (ESS > 5% of total draws). We check for divergent transitions, which would indicate unreliable posterior exploration. Posterior predictive checks compare the distribution of observed data to data simulated from the fitted model, using 100 posterior predictive draws for density overlays (Figure S1) and 1000 draws for computing Bayesian p-values for skewness and kurtosis, computed as P(T(y_rep) ≥ T(y_obs)), where values between 0.05 and 0.95 indicate adequate fit (Gelman et al. 2013: 145). Model comparison between Student-t and Gaussian likelihoods is performed via leave-one-out cross-validation (LOO-CV) using the loo package (Vehtari et al. 2017), with the expected log predictive density (ELPD) as the comparison criterion.

3.6 Inference Targets

Three quantities of interest are derived from the posterior: (1) site mean differences (Δμ = μ_Zhengzhou − μ_Panlongcheng) for each isotope ratio, with 95% credible intervals (CIs); (2) variance ratios (σ²_Panlongcheng / σ²_Zhengzhou) for each ratio, to assess whether the two assemblages exhibit comparable within-site variability; and (3) residual correlations between isotope ratios after accounting for site effects.

4 Results

The raw data distributions for the two assemblages are shown in Figure 3, which displays raincloud plots for each isotope ratio by site. Both assemblages are dominated by highly radiogenic compositions, with a tail extending toward lower, less radiogenic values (mild left-skew) in all three ratios. This asymmetry is most pronounced at Zhengzhou and is slight at Panlongcheng.

4.1 Model Selection

The Student-t model is strongly preferred over the Gaussian alternative. LOO-CV yields ΔELPD = −66.4 (SE = 13.1) comparing Student-t relative to Gaussian (Table 2; the supplementary materials). The negative value indicates superior predictive performance for the Student-t model, and the magnitude exceeds five times its standard error. The model produces zero problematic Pareto k values, whereas the Gaussian model yields four observations with k > 0.7, indicating poor predictive performance for those data points (Figure S2; see Table S4 for detailed LOO-CV results). The estimated degrees-of-freedom parameter median is 1.99 (95% CI: 1.34–3.06; Figure 4), confirming that the lead isotope data exhibit substantially heavier tails than a normal distribution. This finding is consistent with Baxter’s (1999) observation that non-normality is pervasive in lead isotope datasets.

4.2 Convergence

All parameters achieve satisfactory convergence. The maximum \hat{\mathrm{R}} across all parameters is 1.002, well below the conventional threshold of 1.01 (Table S3). The minimum effective sample size ratio is 0.47 (for site-specific dispersion parameters), comfortably exceeding the 5% threshold. Zero divergent transitions occur across all chains, indicating reliable posterior exploration.

4.3 Site Mean Differences

The posterior site means are: Panlongcheng ²⁰⁶Pb/²⁰⁴Pb = 18.96 (95% CI: 18.26–19.63), Zhengzhou ²⁰⁶Pb/²⁰⁴Pb = 20.71 (95% CI: 20.17–21.22); Panlongcheng ²⁰⁷Pb/²⁰⁴Pb = 15.62 (95% CI: 15.52–15.72), Zhengzhou ²⁰⁷Pb/²⁰⁴Pb = 15.85 (95% CI: 15.77–15.92); Panlongcheng ²⁰⁸Pb/²⁰⁴Pb = 39.27 (95% CI: 38.50–39.99), Zhengzhou ²⁰⁸Pb/²⁰⁴Pb = 41.20 (95% CI: 40.62–41.76) (Table 3; Figure 5).

The posterior distributions of site mean differences (Δμ = μ_Zhengzhou − μ_Panlongcheng) reveal statistically significant differences on all three isotope ratios (Table 3; Figure 6). For ²⁰⁶Pb/²⁰⁴Pb, the mean difference is 1.74 (95% CI: 0.88–2.60). For ²⁰⁷Pb/²⁰⁴Pb, the difference is 0.23 (95% CI: 0.10–0.35). For ²⁰⁸Pb/²⁰⁴Pb, the difference is 1.92 (95% CI: 0.99–2.86). In all three cases, the 95% credible intervals exclude zero, providing strong evidence that Zhengzhou bronzes have systematically higher mean isotope ratios than Panlongcheng bronzes.

4.4 Variance Comparison

Despite the significant differences in means, the two sites exhibit similar within-site variability. Variance ratios (Panlongcheng / Zhengzhou) are, ²⁰⁶Pb/²⁰⁴Pb = 1.63 (95% CI: 0.77–3.31), ²⁰⁷Pb/²⁰⁴Pb = 1.57 (95% CI: 0.74–3.18), and ²⁰⁸Pb/²⁰⁴Pb = 1.56 (95% CI: 0.74–3.16). All 95% credible intervals include 1.0, indicating no statistically significant difference in variance between the two assemblages (Figure 7). These residual standard deviations are estimated after accounting for known measurement error via quadrature addition. The posterior probabilities that Panlongcheng has greater variance than Zhengzhou are 0.90, 0.88, and 0.88 for the three ratios, respectively. This suggests that the Panlongcheng data are more scattered than those of Zhengzhou, though the evidence is suggestive rather than conclusive. Notably, this conclusion is sensitive to the choice of likelihood. Under a Gaussian model, the variance ratios are substantially larger (²⁰⁶Pb/²⁰⁴Pb = 2.34 [1.51, 3.65]; ²⁰⁷Pb/²⁰⁴Pb = 2.08 [1.34, 3.22]; ²⁰⁸Pb/²⁰⁴Pb = 1.71 [1.08, 2.70]) and all 95% credible intervals exclude 1.0 (Table S7). This occurs because the Gaussian model, lacking heavy tails, attributes outlying observations to increased variance rather than to the tail behaviour captured by the Student-t likelihood. The Student-t model’s conclusion of comparable variability is therefore the more conservative assessment.

4.5 Residual Correlations

The residual correlations between isotope ratios, estimated jointly across all site–ratio combinations after accounting for site effects, are extremely high, ²⁰⁶Pb/²⁰⁴Pb–²⁰⁷Pb/²⁰⁴Pb = 0.995 (95% CI: 0.991–0.997); ²⁰⁶Pb/²⁰⁴Pb–²⁰⁸Pb/²⁰⁴Pb = 0.992 (95% CI: 0.985–0.996); ²⁰⁷Pb/²⁰⁴Pb–²⁰⁸Pb/²⁰⁴Pb = 0.990 (95% CI: 0.982–0.995). These near-perfect correlations reflect the well-known geochemical coupling between lead isotope ratios, which share a common ²⁰⁴Pb denominator and are produced by related radioactive decay chains. This result underscores the importance of modelling the ratios jointly rather than treating bivariate scatter plots independently.

4.6 Posterior Predictive Checks

Posterior predictive checks indicate adequate model fit for central tendency and asymmetry (Figure S1). Bayesian p-values for skewness range from 0.49 to 0.57 (Table S5) across all ratio–site combinations, indicating no systematic skewness misfit. Bayesian p-values for kurtosis are at or near 1.0 (0.998–1.000), indicating that the fitted model predicts heavier tails than are observed in the data. This is expected behaviour when the estimated ν ≈ 2 is applied to data that, while non-normal, are not as extreme as a Student-t distribution with two degrees of freedom. Crucially, this kurtosis mismatch does not affect the site mean comparison, as confirmed by the stable mean estimates across the Student-t and Gaussian models.

5 Discussion

5.1 Archaeological Interpretation

The Bayesian multivariate analysis reveals a pattern that is not apparent from visual inspection of scatter plots (cf. Figure 2; Figure 8). Panlongcheng and Zhengzhou Shang City bronzes have statistically significantly different mean lead isotope ratios across all three measured ratios. Zhengzhou bronzes exhibit systematically higher mean values, a difference that, while modest in absolute terms (e.g., Δ²⁰⁶Pb/²⁰⁴Pb = 1.74), is nonetheless clearly resolved by the model (95% CI: 0.88–2.60). The comprehensive summary of all model results is presented in Figure 8.

This finding provides quantitative support for the recycling hypothesis proposed by Liu, S. et al. (2023). Their analysis of metallurgical remains from the Panlongcheng Xiaozui workshop identified a group of less radiogenic lead (²⁰⁶Pb/²⁰⁴Pb ≈ 16.5) that they attributed to recycled Lower Erligang bronzes originally transported south from the Central Plains, probably in the earlier Lower Erligang period. The significant downward shift in Panlongcheng’s mean isotope ratios relative to Zhengzhou is consistent with the admixture of this recycled older material, which carries lower radiogenic signatures, into the Panlongcheng bronze assemblage. The “intermediate” lead isotope group (²⁰⁶Pb/²⁰⁴Pb ≈ 16.9–17.5) that Liu, S. et al. (2023) interpreted as a mixture of low-ratio and highly radiogenic lead components would mechanistically produce the lower assemblage mean that our model detects.

At Zhengzhou, the Upper Erligang assemblage appears to rely more heavily on newly procured highly radiogenic lead, consistent with the site’s role as the capital and its direct access to contemporaneous metal supply networks. The lower prevalence of recycled Lower Erligang material at Zhengzhou during this period may reflect the capital’s preferential access to fresh metal supplies.

The significance of this mean difference must be understood within the broader context of the ongoing debate about highly radiogenic lead (HRL) in Shang bronzes (see Liu, S. et al. 2025; Liu, R. et al. 2018; Jin et al. 2017). As Liu, S. et al. (2025) has recently reviewed, HRL has been the most extensively discussed topic in Chinese bronze provenance research over the past four decades. Over 1,000 lead isotope analyses of Shang bronzes have been published, with the majority categorised as HRL (²⁰⁶Pb/²⁰⁴Pb > 19; Jin et al. 2017; Liu, S. et al. 2018; Ma and Cui 2024). A persistent tendency in the literature has been to treat all HRL as a single, undifferentiated category, a monolithic isotopic signature assumed to derive from a common source controlled by the Shang court (Jin et al. 2017). This assumption has, in practice, encouraged scholars to emphasise the similarities between HRL-bearing assemblages across sites and to downplay or overlook inter-site differences such as those between Panlongcheng and Zhengzhou.

However, this monolithic view has been increasingly challenged. Jin Zhengyao first noted that HRL data from the Sanxingdui bronzes distribute along different trend lines from those of the Central Plains (Jin 2008: 106). Liu, R. et al. (2018) subsequently proposed that HRL may have had multiple origins, used by Shang period cultures across different areas. Wang et al. (2021) as well as Ma and Cui (2024), with a substantially larger corpus of data, have confirmed that the chronological change of HRL data trend lines between the Early and Middle Shang is real.

They proposed that the slope of the linear array in lead isotope plots could serve as a criterion for distinguishing different types of HRL. While most previous discussions assumed that these different trend lines represent isotopically distinct HRL resources, Liu, S. et al. (2025) raised the possibility that they are an artefact of different mixing regimes. These are all valuable insights, but they remain fundamentally qualitative. The assessment of whether two slopes or intercepts differ relies on visual judgement rather than on a formal statistical test. It calls for methods capable of detecting subtle distributional differences within the broad HRL category.

The Bayesian multivariate framework introduced in this paper provides precisely such a quantitative test. By modelling the joint distribution of all three isotope ratios simultaneously and estimating site-specific means and variances, the model captures differences in both the position and the dispersion of assemblages along the shared HRL trend. The significant mean differences detected between Panlongcheng and Zhengzhou (Table 3; Figure 6), with Zhengzhou exhibiting higher ²⁰⁶Pb/²⁰⁴Pb, ²⁰⁷Pb/²⁰⁴Pb, and ²⁰⁸Pb/²⁰⁴Pb, formalise what visual approaches have begun to suggest that not all HRL-bearing assemblages are isotopically equivalent, even when they fall along the same broad linear trend. Crucially, our approach quantifies the magnitude of these differences with full posterior uncertainty, replacing subjective assessments with reproducible probabilistic inference.

The comparison in variances between the two sites (Figure 7) is also informative. Both assemblages exhibit comparable within-site scatter, suggesting that the mixing complexity, the number and diversity of metal sources contributing to each site’s bronze production, was broadly similar. This is consistent with the view that Panlongcheng operated a workshop of comparable sophistication to those at the capital, even if the specific metal inputs differed. The slightly higher variance of Panlongcheng may even suggest that its metal sources were more diverse than those of Zhengzhou. Pure reliance on recycled imported bronzes would theoretically produce lower variance, as the recycling process homogenises isotopic signatures.

Taken together, these results support the hypothesis that Panlongcheng had access to partially different metal sources than Zhengzhou. This conclusion aligns with the trace element evidence (Liu, R. et al. 2019), the significantly higher tin content at Panlongcheng (Liu, R. et al. 2017), the distinct crucible technology (Liu, S. et al. 2020), and the evidence for recycling of older bronzes (Liu, S. et al. 2023). The picture that emerges is not one of complete independence but of a regional centre that, while participating in a shared metal circulation network, also drew on local and recycled metal resources.

5.2 Reconciling with Previous Visual Comparisons

Previous visual comparisons of lead isotope scatter plots have generally emphasised the similarities between Panlongcheng and Zhengzhou, noting that both sites yield bronzes with highly radiogenic lead that plot along a common linear trend (Figure 2). This emphasis on similarity is itself a product of the monolithic treatment of HRL discussed above. When all highly radiogenic data are assumed to represent the features of original metal source, inter-site differences are easily attributed to noise rather than to genuine signal of different metal mixing practices. The extremely high residual correlations estimated by our model (≈ 0.99 between all ratio pairs; the supplementary materials) explain why this impression arises. The two sites lie along a shared line in isotope space, and the visual dominance of this linear structure obscures the shift in mean position along the line.

This is a key insight. The difference between the two sites is not one of distinct isotopic clusters but of different average positions along a shared trend (Figures 3 and 8). Visual inspection of scatter plots is poorly suited to detecting such shifts, particularly when within-site scatter is substantial and sample sizes are moderate. The Bayesian model, by marginalising over the nuisance correlations and estimating site means with full uncertainty, resolves a signal that is present in the data but invisible to the eye.

5.3 Applying Bayesian regression to lead isotope data

This paper introduces a Bayesian multivariate regression framework that addresses several long-standing methodological gaps in archaeological lead isotope research. First, by modelling all three isotope ratios simultaneously with estimated residual correlations, it captures the joint distribution of the data rather than treating bivariate scatter plots independently. Second, by using a Student-t likelihood with estimated degrees of freedom, it accommodates the heavy tails that Baxter (1999) showed are characteristic of lead isotope data. Third, by incorporating known analytical measurement error, it responds to Crema’s (2025) call for measurement-error models in non-chronological archaeological applications. Fourth, by producing full posterior distributions for all quantities of interest, it provides probabilistic answers with credible intervals and posterior probabilities rather than binary significance decisions.

The approach offers distinct advantages over both visual biplot inspection and standard frequentist tests. Unlike scatter plots, it is reproducible and quantitative. Unlike classical tests such as Hotelling’s T² or MANOVA, it does not assume multivariate normality, incorporates measurement error, and provides full posterior distributions rather than p-values. The framework is readily transferable to other assemblage comparisons whether between sites, regions, or chronological periods, and can be extended to include additional covariates or hierarchical structure as the research question demands.

More broadly, this study represents a development in the growing application of Bayesian regression methods to compositional and isotopic data in archaeology. While Bayesian approaches have become established in archaeological inference generally (Otárola-Castillo et al. 2022), including radiocarbon chronology (Bronk Ramsey 2009; Banks et al. 2019; Rotunno and Crema 2025), compositional analysis (Vieri et al. 2025), zooarchaeological assemblage modelling (Fernée and Trimmis 2021; Gerbault et al. 2016; Wolfhagen 2024; Ragno 2024), and strontium isotope mobility studies (Scaffidi and Knudson 2020), their adoption for lead isotope provenance research has lagged behind. Our approach addresses all three challenges of this application (strong correlation, heavy tail and measurement error) simultaneously within a single coherent framework.

We note that this approach complements rather than replaces existing methods for lead isotope interpretation. Source attribution methods such as kernel density estimation (De Ceuster et al. 2023), Bayesian mixing models (Sun, Z. et al. 2023b), and multivariate clustering (Tomczyk and Żabiński 2023; Wood and Liu 2023) address the different question of which geological sources contributed to an assemblage.

5.4 Future Perspectives

Several limitations should be noted. First, the known measurement error is specified at 0.1% relative standard deviation, a conservative estimate that may not perfectly reflect the analytical uncertainty of all measurements in the compiled dataset.

Second, sensitivity analyses indicate that the site mean estimates are robust both to the choice of prior on the means (Figure S3) and to the choice of likelihood (Student-t vs. Gaussian; Table S6), whereas the dispersion (sigma) parameters are more sensitive to the likelihood (up to 34% shift; Table S6), reflecting the strong influence of the heavy-tailed likelihood on variance estimation.

Third, the kurtosis posterior predictive check yields p-values of 1.0 across all subgroups, indicating that the model predicts heavier tails than are observed in the data (Figure S1). This mismatch arises because the Student-t distribution is a continuous statistical approximation of an inherently discrete physical process. Bronze compositions result from mixing specific solid end-members, highly radiogenic fresh lead and recycled older bronzes, in a crucible, producing distributions that are bounded by the isotopic compositions of the source materials. The continuous Student-t distribution with ν ≈ 2 predicts outliers beyond these physical bounds, explaining why the model over-estimates tail weight while still capturing the heavy central tendency far better than a Gaussian. This mismatch affects the tails of the predictive distribution but not the site mean comparison, which is well-validated by LOO-CV.

The sample sizes are moderate (38 and 33), and larger assemblages would yield tighter credible intervals and more robust conclusions. In a sense, this demonstrates that extensive analysis of a single assemblage is important and archaeologically meaningful. Previously, a common assumption in lead isotope studies is that a few samples that can “represent” the general distribution range of an assemblage would be enough for a project. However, the approach suggested by the current research shows that 3–5 samples are far from being able to capture the complete distribution information (c. mean, variance, correlation, tail) and a much larger sample size is necessary if a fully quantitative comparison between assemblages is expected.

More broadly, this framework encourages a shift from individual artefact-to-source attribution toward assemblage-level distributional comparison. Traditional approaches that assign individual artefacts to geological sources provide a static view of metal provenance, whereas the rich information contained in the distribution patterns of assemblage data, means, variances, correlations, and tail behaviour, offers a more nuanced understanding of metal circulation dynamics. Future studies should explore hierarchical extensions of this framework to compare multiple sites simultaneously and to incorporate temporal covariates.

6 Conclusions

This study introduces a Bayesian multivariate Student-t regression framework for formally comparing lead isotope ratio distributions between archaeological assemblages. Applied to 71 bronzes from Panlongcheng and Zhengzhou Shang City dating to the Upper Erligang period, the model reveals statistically significant differences in mean lead isotope ratios on all three measured ratios (²⁰⁶Pb/²⁰⁴Pb, ²⁰⁷Pb/²⁰⁴Pb, ²⁰⁸Pb/²⁰⁴Pb), with 95% credible intervals excluding zero in each case. Zhengzhou bronzes exhibit systematically higher mean ratios than Panlongcheng bronzes.

These results provide quantitative support for the hypothesis that Panlongcheng accessed partially different metal sources than the Erligang capital. The lower mean isotope ratios at Panlongcheng are consistent with the admixture of recycled Lower Erligang bronzes. The similar variances between the two sites suggest comparable mixing complexity at both workshops, while the extremely high residual correlations between ratios explain why visual scatter plot comparisons have emphasised similarity rather than difference.

Methodologically, this framework addresses several recognised limitations of current practice in archaeological lead isotope studies. It replaces subjective visual assessment with reproducible probabilistic inference, and models all three isotope ratios jointly. This method also accommodates heavy-tailed distributions, and it incorporates known analytical measurement error. The approach is readily transferable to other assemblage comparisons and provides a template for the quantitative comparative analysis that the growing lead isotope database increasingly demands.

References

[1] Albarède Desaulty , F. Blichert-Toft . A Geological Perspective on the Use of Pb Isotopes in Archaeometry. Archaeometry, 2012, 54(5): 853–867

[2] Banks W. E, P. Bertran, S. Ducasse, L. Klaric, P. Lanos, C. Renard , M. Mesa. . An Application of Hierarchical Bayesian Modeling to Better Constrain the Chronologies of Upper Paleolithic Archaeological Cultures in France between ca. 32,000–21,000 Calibrated Years before Present. Quaternary Science Reviews, 2019, 220: 188–214

[3] Baxter M. J. . On the Multivariate Normality of Data Arising from Lead Isotope Fields. Journal of Archaeological Science, 1999, 26(1): 117–124

[4] Bronk Ramsey C. . Bayesian Analysis of Radiocarbon Dates. Radiocarbon, 2009, 51(1): 337–360

[5] Bürkner . brms: An R Package for Bayesian Multilevel Models Using Stan. Journal of Statistical Software, 2017, 80: 1–28

[6] Cao D., E. R. Crema , E. Pomeroy. . Estimating Intralimb Proportions for Commingled Remains. International Journal of Osteoarchaeology, 2024, 34(5): e3326

[7] Chen J., J. Zhang, Q. Fang, C. Zhang, C. Gao , S. Chen. . Jin Dao Xi Hang: Metal Circulation Network in the Middle Yangtze River Region during the Shang and Zhou Periods金道锡行——简论商周时期长江中游地区金属流通网络. Jianghan Kaogu, 2024, 2024(5): 11–23

[8] Crema E. R. . Statistical Modelling in Archaeology: Some Recent Trends and Future Perspectives. Journal of Archaeological Science, 2025, 180: 106295

[9] De Ceuster S., D. Machaira , P. Degryse. . Lead Isotope Analysis for Provenancing Ancient Materials: A Comparison of Approaches. RSC Advances, 2023, 13: 19595–19606

[10] Fernée L. , C. P. Trimmis . Detecting Variability: A Study in Applying Bayesian Multilevel Modelling to Archaeological Data. Journal of Archaeological Science, 2021, 128: 105346

[11] Gelman, A., J. B. Carlin, H. S. Stern, D. B. Dunson, A. Vehtari, and D. B. Rubin. 2013. Bayesian Data Analysis, 3rd ed. Boca Raton: CRC Press.

[12] Gerbault P., R. Gillis, J.-D. Vigne, A. Tresset, S. Bréhard , M. G. Thomas. . Statistically Robust Representation and Comparison of Mortality Profiles in Archaeozoology. Journal of Archaeological Science, 2016, 71: 24–32

[13] Hubei Provincial Institute of Cultural Relics and Archaeology. 2001. The Panlongcheng site: Report of Archaeological Excavation from 1963-1994盘龙城——1963~1994年考古发掘报告. Beijing: Wenwu Chubanshe (in Chinese).

[14] Jin, Z. 2008. Lead Isotope Archaeology In China中国铅同位素考古. Hefei: University of Science and Technology of China Chubanshe (in Chinese).

[15] Jin Z., R. Liu, J. Rawson , A. M. Pollard. . Revisiting Lead Isotope Data in Shang and Western Zhou Bronzes. Antiquity, 2017, 91(360): 1574–1587

[16] Lewandowski D., D. Kurowicka , H. Joe. . Generating Random Correlation Matrices Based on Vines and Extended Onion Method. Journal of Multivariate Analysis, 2009, 100(9): 1989–2001

[17] Li, H., 2020. Resources and Society: Metal Circulation in the Shang and Zhou Periods 资源与社会:以商周时期铜器流通为中心. Beijing: China Social Sciences Press (in Chinese).

[18] Liu, L., and X. Chen. 2003. State Formation in Early China. London: Duckworth.

[19] Liu, L., and X. Chen. 2012. The Archaeology of China: From the Late Paleolithic to the Early Bronze Age. Cambridge: Cambridge University Press.

[20] Liu R. , A. M. Pollard. . Asking Different Questions: Highly Radiogenic Lead, Mixing and Recycling of Metal and Social Status in the Chinese Bronze Age. Mineralogical Magazine, 2022, 86(4): 677–687

[21] Liu R., A. M. Pollard, J. Rawson, X. Tang, P. Bray , C. Zhang. . Panlongcheng, Zhengzhou and the Movement of Metal in Early Bronze Age China. Journal of World Prehistory, 2019, 32(4): 393–428

[22] Liu R., A. M. Pollard, J. Rawson, X. Tang , C. Zhang. . . Jianghan Kaogu, 2017, ((3)): 119–129

[23] Liu R., J. Rawson , A. M. Pollard. . Beyond Ritual Bronzes: Identifying Multiple Sources of Highly Radiogenic Lead across Chinese History. Scientific Reports, 2018, 8: 11770

[24] Liu S., Z. Sun, J. Zhang, J. Lin, K. Chen, J. Chen, J. Cui , J. Mei. . 50 Years of Bronze Provenance Studies: A Perspective from China. Journal of Archaeological Science, 2025, 180: 106272

[25] Liu S., K. Chen, T. Rehren, J. Mei, J. Chen, Y. Liu , D. Killick. . Did China import metals from Africa in the Bronze Age?. Archaeometry, 2018, 60(1): 105–117

[26] Liu S., Q. Zou, J. Lu, K. Chen , J. Chen. . The Provenance of Bronzes from the Site of Xiaozui at Panlongcheng盘龙城遗址小嘴金属物料溯源研究. Jianghan Kaogu, 2023, (4): 131–138

[27] Liu S., Q. Zou, J. Lu, K. Chen , J. Chen. . Scientific Analysis and Investigation of Metallurgy-related Remains of the Shang Dynasty at the Xiaozui of the Panlongcheng Site盘龙城遗址小嘴商代冶金遗物的分析与研究. Jianghan Kaogu, 2020, (4): 126–137

[28] Ma R. , J. Cui. . Shang and Beyond: Preliminary Exploration on the Evolution and Circulation Pattern of the Highly Radiogenic Lead Resources in the Shang Dynasty商与远方——商代高放射性成因铅资源时空演变与流通模式的初步探索. Journal of National Museum of China, 2024, (9): 84–104

[29] Nan P., Y. Qin, T. Li , Y. Dong. . Analysis of the Casting Origins of Some Shang Dynasty Bronzes Unearthed at Panlongcheng, Hubei湖北盘龙城出土部分商代青铜器铸造地的分析. Wenwu, 2008, (8): 77–82

[30] Otárola-Castillo G. Torquato , E.R. Wolfhagen , M. Hill Jr. , J. Buck . Beyond Chronology: Using Bayesian Inference to Evaluate Hypotheses in Archaeology. Advances in Archaeological Practice, 2022, 10(4): 397–413

[31] Papageorgiou I. . Ceramic Investigation: How to Perform Statistical Analyses. Archaeological and Anthropological Sciences, 2020, 12: 210

[32] Peng, Z., Z. Wang, W. Sun, S. Liu, and X. Chen. 2001. “Lead Isotope Analysis of Shang Dynasty Bronze Artefacts from Panlongcheng 盘龙城商代青铜器铅同位素示踪研究.” In The Panlongcheng site: Report of Archaeological Excavation from 1963-1994盘龙城:1963年—1994年考古发掘报告. Beijing: Wenwu Chubanshe, pp. 552–558 (in Chinese).

[33] Ragno R. . Sheep and Goats Taxonomic Abundance Trends in 1st Millennium CE Southern Italy: Multilevel Bayesian Modelling of NISP Datasets. Journal of Archaeological Science, 2024, 171: 106068

[34] Rotunno R. , E. R. Crema. . Bayesian Analyses of Radiocarbon Dates Suggest Multiple Origins of Ceramic Technology in Early Holocene Africa. Nature Communications, 2025, 16: 8819

[35] Scaffidi B. K. , K. J. Knudson. . An Archaeological Strontium Isoscape for the Prehistoric Andes: Understanding Population Mobility through a Geostatistical Meta-analysis. Journal of Archaeological Science, 2020, 117: 105121

[36] School of History of Wuhan University , Hubei Provincial Institute of Cultural Relics , Archaeology Site Museum . Excavation Report on the Xiaozui Locus at Panlongcheng, Wuhan, 2015–2017武汉市盘龙城遗址小嘴 2015–2017年发掘简报. Jianghan Kaogu, 2019, (6): 15–34

[37] Stan Development Team. 2023. RStan: The R Interface to Stan. R package version 2.32.

[38] Sun, S., R. Han, T. Chen, T. Saito, M. Sakamoto, and I. Taguchi. 2001. “Lead Isotope Analysis of Bronze Artefacts from Panlongcheng 盘龙城出土青铜器的铅同位素比值测定报告.” In The Panlongcheng site: Report of Archaeological Excavation from 1963-1994盘龙城:1963年—1994年考古发掘报告. Beijing: Wenwu Chubanshe, pp. 545–551 (in Chinese).

[39] Sun Z., S. Liu, S. Yang, K. Chen , J. Chen. . Investigating the Origins of Metals Used in the Early Shang Capital of Zhengzhou. Journal of Archaeological Science: Reports, 2023a, 48: 103872

[40] Sun Z., S. Liu, J. Zhang, K. Chen , B. Kaufman. . Resolving the Complex Mixing History of Ancient Chinese Bronzes by Manifold Learning and a Bayesian Mixing Model. Journal of Archaeological Science, 2023b, 151: 105728

[41] Tian, J. 2013. “A Study of Erligang-Period Bronzes Excavated in the Zhengzhou Area 郑州地区出土二里冈期铜器研究.” PhD diss., University of Science and Technology of China中国科学技术大学, Hefei (in Chinese).

[42] Tomczyk C. , G. Żabiński. . A PCA-AHC Approach to Provenance Studies of Non-Ferrous Metals with Combined Pb Isotope and Chemistry Data. Journal of Archaeological Method and Theory, 2023, 31: 93–143

[43] Vehtari A., A. Gelman , J. Gabry. . Practical Bayesian Model Evaluation Using Leave-one-out Cross-validation and WAIC. Statistics and Computing, 2017, 27: 1413–1432

[44] Vieri J., E. R. Crema, M. A. Uribe Villegas, J. Sáenz Samper , M. Martinón-Torres. . Beyond Baselines of Performance: Beta Regression Models of Compositional Variability in Craft Production Studies. Journal of Archaeological Science, 2025, 173: 106106

[45] Wang Q., J. Guo, J. Chen, S. Liu, Z. Fang, M. Li , H. Fang. . Lead Isotope Analysis of Shang Dynasty Bronzes Unearthed at the Liujiazhuang Site, Jinan 济南市刘家庄遗址出土商代青铜器的铅同位素分析. Kaogu, 2021, (7): 106–120

[46] Wolfhagen J. L. . Estimating the Ontogenetic Age and Sex Composition of Faunal Assemblages with Bayesian Multilevel Mixture models. Journal of Archaeological Method and Theory, 2024, 31: 507–556

[47] Wood J. R. , Y. Liu. . A Multivariate Approach to Investigate Metallurgical Technology: The Case of the Chinese Ritual Bronzes. Journal of Archaeological Method and Theory, 2023, 30: 707–756

[48] Zhang, C. 2014. “Erligang: A Perspective from Panlongcheng.” In Art and Archaeology of the Erligang Civilization, edited by K. Steinke and D. C. Y. Ching. Princeton University Press, pp. 51–63.

[49] Zhang C. . On Bronzes of the Erligang Culture Period: General Characteristics and Significance论二里冈文化时期青铜器:一般特性及意义. Jianghan Kaogu, 2024, (3): 83–88

Rights & permissions

The Author(s) 2026. This article is available under open access at journal.hep.com.cn.