成人骨骼年龄估计的完全贝叶斯方法：基于马尔可夫链蒙特卡洛采样的多变量潜在特质建模

《American Journal of Biological Anthropology》：A Fully Bayesian Approach to Adult Skeletal Age Estimation: Multivariate Latent Trait Modeling With Markov Chain Monte Carlo Sampling

　　

摘要

目的

本文考察贝叶斯转变分析在改进骨骼学数据集死亡年龄估计中的潜力，并采用 NIMBLE 与 JAGS（Just Another Gibbs Sampler）软件实现。研究人员旨在探索该方法对于不同骨骼数据集及不同年龄标志物的适用性。

材料与方法

转变分析基于转变年龄（age-at-transition）与总体先验（population priors）对骨骼年龄标志物进行建模估计。本研究通过将多个潜在特质（latent traits）纳入贝叶斯框架，对该方法学进行了扩展。所建 probit 回归模型可直接估计死亡年龄（age-at-death）与 Gompertz 分布参数，且无需依赖预先设定的转变年龄参考总体或总体先验。分析纳入了来自 3 个不同来源的数据，这些数据包含不同组合的年龄标志物，且其死亡年龄均已知。

结果

该模型在各数据集上均表现良好，并成功估计了死亡模式（mortality patterns）。样本量确有影响，但即使较小样本仍可产生有用结果；特质数量亦然。简单的多重 probit 回归会低估年龄区间，但经校准后性能得到改善。即使存在缺失数据，模型仍可产生稳健估计。

讨论

结果表明，该模型能够在不同总体中有效改进死亡年龄估计及死亡模式推断。该贝叶斯转变分析为生物考古学（bioarchaeology）应用中的年龄估计提供了一种可靠且灵活的方法，尤其适用于形态学年龄指征质量不一的数据集。总体间的相容性甚至不是必需条件。唯一必要的是：对每个骨骼组合（skeletal assemblage）具有观察者内部一致的年龄标志物评定，并且所讨论特质必须随生物学年龄单调变化。

该文发表于《American Journal of Biological Anthropology》，围绕成人骨骼死亡年龄估计这一长期难题，提出了一种尽可能符合“Rostock Manifesto”要求的完全贝叶斯（Bayesian）分析框架。研究背景在于，传统骨骼年龄估计通常依赖已知死亡年龄的参考样本，将某些骨形态阶段直接对应到特定年龄区间；然而，这类方法在古人口学（paleodemography）中存在根本局限，因为研究真正需要的是在已知骨骼形态证据条件下个体死于某年龄的后验概率，而这一概率又依赖目标总体死亡年龄分布。换言之，若不先理解目标总体的死亡模式（mortality pattern），单纯套用参考人群的年龄分期便可能产生系统偏差。既往许多方法虽冠以“贝叶斯”之名，但往往仍将年龄转变阈值、总体先验或参考分布视为固定已知量，因此并未真正将所有未知参数作为随机变量处理。正是在这一问题意识下，研究人员尝试构建一种不依赖外部参考总体、可在单一模型中同时估计年龄转变、总体死亡分布及个体死亡年龄的新方案。

研究人员开展的核心工作，是把传统转变分析中分离的若干步骤整合进一个统一的贝叶斯潜在特质（latent trait）模型中，并通过马尔可夫链蒙特卡洛（Markov Chain Monte Carlo，MCMC）采样联合估计参数。文章以有序 Probit 回归（ordered probit regression）为基础，将骨骼年龄指征的离散等级视为潜在连续过程的外在表现；每一性状等级之间的阈值并非预先给定，而是在模型中一并估计。与此同时，个体死亡年龄服从 Gompertz 分布，并直接由目标样本数据推断其参数。该设计的意义在于，只要相关骨学性状随生物学年龄单调变化，就不必再要求外部参考样本提供固定的转变年龄或总体年龄先验，从而减少跨总体不相容、年龄模仿（age mimicry）及观察标准迁移带来的误差。论文最终表明，这一方法在单一性状、多性状及存在缺失值的不同场景下均具有可用性，并能同时提供个体年龄估计与总体死亡模式推断，因此对于生物考古学中保存状况不一、指标来源复杂的骨骼组合具有重要方法学意义。

方法上，研究人员主要使用 3 类关键技术路径。其一，构建单变量或多变量有序 Probit 潜在特质模型，并分别采用条件独立的 multiple-normal 模型与考虑性状相关性的 multivariate normal 模型。其二，以 Gompertz 分布描述成人死亡年龄分布，在模型内同步估计死亡风险参数。其三，采用 R、JAGS 与 NIMBLE 实施 MCMC 采样，并以潜在尺度缩减因子（PSRF）和有效样本量（ESS）评估收敛质量。样本来源包括：Konigsberg 等提供的 “Allmales” 已知死亡年龄男性耻骨联合数据，Navega 等公开的 Coimbra Identified Skeletal Collection（CISC）与 XXI-ISC 组成的 CAMSAD 多性状数据，以及 Christ Church, Spitalfields 的髂骨耳状面数据；另通过抽样、增加性状数与人为引入缺失值测试模型稳健性。

5.1 “Allmales” Dataset

在“Allmales”数据集中，研究人员以 Suchey–Brooks 耻骨联合 6 阶段评分为基础，检验不同样本量对模型表现的影响。结果显示，当样本量为 200 时，模型得到的平均年龄略高于真实死亡年龄，偏差（Bias）为 2.96 岁；残差年龄斜率（residual age slope）为 0.21，说明年轻个体倾向被高估、年长个体倾向被低估，但已知年龄与估计年龄之间仍具有较高相关性。年龄估计的平均绝对误差约为 9.88 岁，95% 最高密度区间（highest density interval，HDI）中位宽度约为 50.72 岁，覆盖率为 90%，略低于理论 95%。研究人员进一步比较了已知死亡年龄计算出的阶段转变阈值与模型估计值，发现二者仅相差数月到数年，且各阈值的 95% 可信区间均覆盖经验值，表明模型对阶段转变过程具有较好恢复能力。

在缩小样本量后，100、50 与 25 例样本的各项拟合指标并未出现显著恶化，提示该方法并不一定要求极大样本才能产生有意义结果；只有当样本量降至 10 时，准确性与稳定性才明显下降。研究人员还构造了一个更接近 Gompertz 总体的低死亡率样本进行比较。该样本的估计准确性与原始样本相近，但覆盖率更接近理论值，说明若目标人群的死亡结构更符合 Gompertz 分布，则模型的区间表现会更理想。文章同时指出，若使用均值而非众数（mode）作为个体年龄点估计，偏差与残差年龄斜率会发生明显变化，提示年龄估计质量评价高度依赖点估计表达方式。

5.2 CAMSAD

在 CAMSAD 数据集中，研究人员考察了随着年龄性状数目增加，模型性能如何变化。该数据集涵盖多个解剖区域的 64 个性状，是三组材料中形态信息最丰富的一组。由于完整多变量模型在此数据规模下计算代价过高，研究采用“简单”多重有序 Probit 回归，并在最终结果中加入 Gaussian 噪声进行不确定性校准。结果显示，随着性状由 5 个逐步增至 64 个，Bias、残差年龄斜率、平均绝对误差、均方根误差（RMSE）和连续排序概率评分（CRPS）均明显改善；尤其在 64 性状条件下，Bias 接近 0，残差年龄斜率仅约 0.02–0.04，平均绝对误差约 6.4–6.5 岁，表明多性状整合显著提升了年龄估计性能。

然而，研究同时发现，若简单假定各性状条件独立，则尽管 95% HDI 范围可压缩至约 10 年，但真实覆盖率仅约 55%，明显过于乐观，说明忽略性状间相关性会严重低估不确定性。经加入标准差 7.5 岁的校准噪声后，95% HDI 扩大至约 30.68 岁，覆盖率提高到 91.5%，且多数覆盖水平的精确二项检验不再显著。这一结果证明，多重 Probit 回归虽然能够利用大量年龄相关信息提高点估计性能，但若不对相关结构或后验不确定性进行处理，其区间估计会失真。由此，文章强调校准对于高维但条件独立假设不成立的数据尤为重要。

5.3 Spitalfields

在 Spitalfields 数据集中，研究人员使用 Buckberry and Chamberlain 髂骨耳状面方法的 5 个性状，采用 multivariate normal 模型直接处理性状间相关性。结果显示，该模型对总体死亡模式拟合较好，Gompertz 参数能够较准确再现真实年龄分布。点估计方面，Bias 为 3.18 岁，残差年龄斜率为 0.30，平均绝对误差为 11.68 岁，RMSE 为 14.92 岁；95% HDI 的中位宽度为 48.42 岁，覆盖率为 92.22%，与理论值较为接近。与“Allmales”类似，该结果也表明基于单一骨学指标体系的成人年龄估计仍保有较宽的不确定区间，但其总体覆盖表现是可接受的。

研究人员随后通过删除 20%、30%、40% 和 50% 的性状记录，模拟考古骨骼材料常见的缺失情境，并改用 simple multiple model 处理。结果表明，20% 缺失对拟合质量影响较小，30% 和 40% 缺失时准确性下降但仍在可接受范围，而 50% 缺失时偏差和误差明显加大。进一步分析指出，某些情况下低等级阶段被完全删除，会导致模型对年龄结构的推断系统性偏移。文章还比较了耳状面各性状阈值的经验转变年龄与模型估计值，发现大多数阈值较为接近，但个别性状如 MA 的某一阈值出现异常，且对应参数收敛质量不足，提示特定性状的分布结构可能限制模型表现。

6 Discussion

讨论部分围绕准确性、系统偏差、年龄区间与既往研究比较展开。首先，在准确性层面，研究显示样本量更大、性状更多、缺失更少通常会提升表现，但不同性状体系的收益并不均一。CAMSAD 多性状系统最佳，Suchey–Brooks 次之，Buckberry–Chamberlain 在本研究情境下略逊。文章特别指出，CRPS 作为基于完整后验分布的全局指标，比仅考察观测点概率密度的 TMNLP 更能稳定反映样本型年龄估计模型的性能。

其次，在 Bias 与残差年龄斜率方面，最佳情形下 Bias 可控制在 ±3 岁以内，甚至接近 1 岁，但年轻被高估、老年被低估的问题并未被完全消除。研究认为，较多性状和较充分的阶段覆盖有助于减轻这一“回归向均值”现象，但更多性状并不自动保证更低的残差年龄斜率，性状本身与年龄的相关强度及分期结构同样关键。

再次，在年龄区间、精度与效率方面，文章强调年龄估计本质上是“区间宽度”与“真实覆盖率”的权衡。Suchey–Brooks 与 Buckberry–Chamberlain 两类传统指标所导出的可信区间通常接近或超过 50 年，而 CAMSAD 在充分校准后约为 30 年，说明多性状整合确可提高区间效率，但前提是要合理处理相关性和不确定性。与其他已发表研究相比，本文模型在多数组合下至少达到相当水平，某些指标如 CAMSAD 的残差年龄斜率甚至更优，但作者也强调仍需更多验证研究支持。

文章还系统归纳了该方法的优势与局限。优势在于：无需“恒定性假设”（invariance hypothesis），不依赖参考总体转变年龄与总体先验，可避免年龄模仿，并能同步估计 Gompertz 死亡风险参数；即使样本仅 25 例，也可能得到有用结果。局限在于：模型默认成人死亡年龄服从 Gompertz 分布，因此不适于儿童全年龄段；灾变性死亡样本或近均匀年龄结构样本可能削弱拟合；多变量有序回归难以直接处理缺失值；若性状之间高度相关或与年龄相关性过弱，也会降低模型效果；此外，该方法更适用于骨骼群体研究，而非单一个体的法医学场景。

研究结论部分指出：研究人员认为，遵循“Rostock Manifesto”精神的贝叶斯版转变分析是可行的。该方法虽不能解决成人死亡年龄估计的全部问题，但至少能够缓解其中一部分关键困难。通过把参考总体界定与转变参数确定合并进同一模型，年龄转变参考总体及其先验不再是必需条件；当缺乏合适参考样本时，尤其是在年代较早的考古骨骼组合研究中，这一点具有特殊价值。该方法可应用于任何随年龄增长而变化的骨骼特征，并为对既往资料进行重新分析提供了统一接口。尽管在缺失值处理、相关性分组加权以及随年龄变化的校准机制等方面仍有改进空间，这一方法仍构成当前年龄估计路径之外的一种真正替代方案，尤其适用于保存状态差异显著且可获得多种年龄性状的材料；即便样本量较小，也可能得到有用结果。