祝孔豪 李 斌 王 康 郭钰伦 王维康 徐 军

(1. 中国科学院水生生物研究所,武汉 430072; 2. 内江师范学院,内江 641100; 3. 深圳市利源水务设计咨询有限公司,深圳 528031)

稳定同位素技术是研究食物网生态学的重要技术手段[1—4]。与直接观察摄食或胃肠含物分析等技术手段相比,稳定同位素技术能反映较长时间尺度内消费者营养来源的整合特征[1,5]。基于稳定同位素质量平衡模型,稳定同位素技术可用于消费者营养溯源,即确定多种营养来源对消费者的贡献比重[6—8]。作为稳定同位素质量平衡模型的重要分支,贝叶斯混合模型考虑了营养来源的不确定性,并允许纳入协变量,所以近十多年来得到了广泛的应用[9]。通常使用贝叶斯混合模型来估计不同营养来源的贡献; 此类模型在贝叶斯混合模型输出了每个营养来源对消费者贡献比例的概率分布特征[10],因此其假设之一就是营养来源对消费者贡献比重满足(0,100%)的条件。当引入某种营养来源进入模型时,就存在了消费者确实依赖了这种营养来源的先验假设; 一个营养来源的贡献比例高,必然会降低其他营养来源的贡献比例,导致营养来源贡献比例不确定性增加[11]。因此,评估贝叶斯混合模型性能,就需要评价贝叶斯混合模型拟合结果的优劣,及其与生态学理论和研究背景的匹配程度。

本文基于实测同位素数据集(蒙古鲌Culter mongolicus mongolicus同位素数据集),通过识别消费者营养功能类群特征、改变营养来源先验信息特征,构建系列贝叶斯模型; 通过比较模型总体性能、实测值与预测值差异,及先验信息和后验信息差异等多种模型性能评价方法,来描述模型性能评价的方法和过程,以此为应用稳定性同位素技术开展消费者营养溯源研究,提供模型性能评估体系。

1 稳定同位素数据

1.1 鲌稳定同位素数据集

本研究中蒙古鲌数据引自李斌等[12]。实验材料于2010年7月采自三峡库区腹地北岸支流小江。同位素样品在西南大学地理科学学院地球化学与同位素实验室(设备型号: Flash EA112HE DELT plus XP)完成测试,其分析精度δ13C<0.02‰,δ15N<0.03‰。样品的收集、处理和保存方法如徐军等[13]所描述。

1.2 种群营养功能群

采用K-means聚类对消费者稳定同位素值进行聚类和种群营养功能群划分。K-means聚类属于无监督学习聚类。基本原理是从n个数据对象任意选择k个对象作为初始聚类中心点[14]; 根据每个聚类对象的中心点,计算每个数据对象与这些中心点的距离; 并根据最小距离重新对相应数据对象进行划分; 重新计算每个聚类中心点,反复循环至每个聚类不再发生变化[14,15]。

较为常见的聚类效果评估方法是总体类内误差平方和(Total within sum of the squared errors)评估(公式 1),即肘部法则(Elbow Method)[14,15]。

式中,Ci是第i个簇;p是Ci中的样本点;mi是Ci的质心(Ci中所有样本的均值);SSE是所有样本的聚类误差[14]。随着聚类数k的增大,总体类内误差平方和减小,SSE随着分组增加而降低; 在达到某个临界点时,SSE得到极大改善,之后缓慢下降,形成一个“肘部”形状的关系; 这个临界点附近的聚类数k即可作为聚类性能较好的分组[14,15]。

本研究分析表明,在对消费者蒙古鲌的稳定同位素值进行聚类过程中,发现在k=4、5、6时,SSE值形成的“手肘”形状关系; 但是分组的增加所获得的SSE值下降的回报不高(图 1A)。结合湖泊生态学基本理论和消费者种群内营养功能群的特点[16—18],本文设置了一个聚类k值的选择条件,即某k值(k>1)聚类所获得的SSE值与k=1时SSE值的比值小于15%时,则分组回报不足。基于这个聚类原则,对消费者蒙古鲌的稳定同位素值进行聚类,最佳聚类数选3(图 1A),分别代表鱼类食性为主(Piscivorous)、浮游动物食性为主(Planktivorous)和底栖动物食性为主(Benthivorous)等食物网功能类型(图 1B)。种群营养功能群的识别,既有助于减小消费者同位素变异所带来的不确定,又有助于理解消费者在食物网的功能[19,20]。

图1 种群内营养功能群的确定Fig. 1 Within-population trophic functional groupsA. K-means聚类与总体类内误差平方和百分比特征; B. 鱼类食性为主、浮游动物食性为主和底栖动物食性为主的3个功能类型A. K-means cluster and percentage of total within sum of the squared errors; B. Three trophic functional groups,including piscivorous,planktivorous,and benthivorous

1.3 营养富集因子与源矫正

构建贝叶斯混合模型的一个重要前提条件是消费者与营养来源间的同位素的差异[21,22]。这种差异由消化和代谢过程中同位素分馏引起,被称为营养富集因子(Trophic enrichment factor,TEF,Δ),其定义为Δ =δtissue–δdiet(δ代表样品的同位素比值与标准物质的同位素比值之间的相对差异)[23]。由于TEF受到从生理到营养来源的多种因素影响,实际研究中多采用营养级富集因子的统计平均值(或多个统计平均值)进行分析,以获得最接近富集因子真实值的近似值,从而更好地估算该种类的食源组成[21,22]。例如,δ15N在相邻两个营养级间所产生的富集因子(Δδ15N)在3‰—5‰[24],而δ13C在相邻营养级间富集因子(Δδ13C)为0.4‰—1.0‰[25]。

在本研究中,消费者蒙古鲌是典型的肉食性消费者,本研究使用Δδ15N=(3.4±0.99)‰和Δδ13C=(0.39±1.14)‰对营养来源的稳定碳氮同位素进行矫正,即Δ+δdiet[21,22]。结合蒙古鲌食性对食物来源进行整合,形成四类营养来源用于模型构建,包括浮游动物、底栖动物、浮游鱼类和底栖鱼类(图 2,食物来源δ值数据来自李斌等[12]的实测数据)。TEF的不确定性(标准误差),通过统计平方公差法(Root-Sum-Squares Error),整合到食物来源同位素的误差中(公式2),用来反映数据误差的整体特征[26]。

式中,σ为标准差。

同位素混合空间中(图 2)营养来源和消费者同位素之间的共线性,及不同营养来源之间的差异水平等情况,可会导致统计上多种营养来源贡献的等效解决方案,也需尽量避免[27]。结合多元方差分析(Multivariate analysis of variance,MANOVA),结果表明4种潜在食物来源表现出一种或一种以上的稳定同位素显著差异(P<0.05),且共线性特征不显著(图 2),适合进一步分析消费者的食物来源。

1.4 建模数据质量检验

TEF的选择及源矫正的合理性,对稳定同位素质量平衡混合模型的结果影响很大[21]。混合模型的基本假设之一就是,消费者的稳定同位素必须属于多种食物来源稳定同位素所定义的同位素混合空间[9]。特别需要指出的是,本研究案例中的贝叶斯混合模型,其数学性质是即使数据不符合稳定同位素混合的基本假设,也会获得方程的解[9,10]。因此,必须检查消费者稳定性同位素数据是否绝大多数落入多种营养来源确定的同位素混合空间中[28]。

尽管通过营养富集因子校正后的稳定同位素混合空间(图 2),可以初步进行判断; 但直观观察很难准确确定哪些样本属于混合多边形。由于贝叶斯建模包含不确定性(本例中为均值和标准偏差),需要借助一定的统计学方法来判断哪些消费者在模型求解过程中产生错误的风险较高。本研究采取混合多边形迭代模拟的方法来判别数据落入混合多边形的可能性。基本方法是,基于TEF校正后的各种营养来源的同位素均值和标准偏差,迭代生成10000次稳定同位素混合多边形; 进一步计算消费者稳定同位素落入这些混合多边形的频次; 消费者稳定同位素落在>0.05可能性区域,为数据质量符合建模需要(图 3A、3B和3C)。

图2 营养富集因子校正后的稳定同位素混合空间Fig. 2 Trophic enrichment factor corrected isospace

如前文所述,同位素混合空间中营养来源和消费者之间的共线性,可为营养来源贡献提供多种在统计上等效的解决方案[9,10]。引起这种不确定性的主要特点是,当消费者的稳定同位素处于由营养来源所构成的同位素空间中心区域时,模型无法确定营养来源比例。因此,当稳定同位素混合空间质心区域消费者样本过多时,共线性增加,稳定同位素混合空间模型求解不会收敛或预测值与观测值不能较好匹配[29]。本研究采取混合多边形迭代模拟的方法来判别数据质量。基本方法是,基于TEF校正后的各种营养来源的同位素均值和标准偏差,迭代生成10000次高风险混合同位素空间(以质心为中心的50%不规则多边形面积内)[28,29],进一步计算消费者同位素落入风险高的混合空间的次数; 通过计算落入混合空间次数与迭代次数的比值,基于统计学频次计算落入高风险稳定同位素混合空间的概率,来检验数据建模的质量; 消费者稳定同位素落在<0.95可能性区域,为数据质量符合建模需要(图 3D、3E和3F)[28,29]。

由图 3所示(等值线颜色深浅显示了概率轮廓),本研究数据总体质量较高。图 3A、3B和3C则显示了消费者蒙古鲌同位素值的变化将如何影响营养来源混合模型合理求解的概率[28]。95%概率轮廓内的消费者蒙古鲌样品可用于混合模型使用,也就是说其中一个样本(16号个体)位于营养来源之外,因此模型解释程度低,在后续建模分析中予以剔除。图 3D、3E和3F则显示了消费者蒙古鲌同位素值的变化将如何影响营养来源混合模型低估风险的概率[29]。消费者蒙古鲌样品落入风险区的概率总体低于50%,未出现高于95%概率的样本。

图3 混合多边形迭代模拟Fig. 3 Mixing polygon simulation

需要强调的是,如果在上述统计方法检验过程中,发现较多的样本不适用于建模,则必须考虑研究假设中是否存在1)测量过程的错误、和/或2)消费者的重要食物来源被忽略、和/或3)营养富集因子使用不合理等问题。

2 贝叶斯质量平衡混合模型

2.1 先验信息

“从先验信息中获益(Profiting from prior information)”是贝叶斯建模中的一个常见想法[30,31]。因此,先验信息的确定是稳定性同位素质量平衡混合模型的一个重要条件。为特定研究设计建立合理的信息先验是建立贝叶斯模型的最佳方案。因此,在开展研究前,需要投入大量的时间和精力进行数据收集和分析,了解哪些数据满足了前期的假定条件,提升进一步建模分析的准确性[32]。当观测样本量或代表性有限时,贝叶斯模型将观测值与先验信息相结合,提升模型预测的准确性。在消费者营养溯源研究中,尽管使用胃肠含物来评估营养来源具有一定的局限性[33],但为同位素技术的研究提供了较好的先验信息,可以提高混合模型预测的准确性[34]。此外,环境中潜在食物来源丰度、生物量和消费者摄食行为习性等,均可作为重要的先验信息,以提高混合模型预测的准确性[34]。

本研究展示了3个营养来源贡献的先验信息,包括默认先验(Uninformative)、信息先验(Informative)和高信息先验(Informative with SD; 图 4)。针对4种营养来源,包括浮游动物、底栖动物、浮游鱼类和底栖鱼类,默认先验的营养来源贡献由均值为零、标准偏差为1的正态分布混合确定,并进行中心化对数比转换(Centralized logarithm ratio transformation)[35],从而产生均值为0.25且营养来源边际贡献大的分布特征(SD≈0.2)。信息先验是基于蒙古鲌的食性分析数据,4种来源质量比例均值分别为16.7%、28.3%、20.6%和34.4%,依据上述方法获得营养来源边际贡献水平较大的分布特征(SD≈0.2)。高信息先验的情景则是在此基础上,结合食性数据变化的标准偏差特点,在均值不变的条件下,对所有营养来源使用了SD ≈ 0.06的标准偏差(图 4)。在进一步建模分析中,本研究采用高信息先验的情景来分析蒙古鲌数据集。

图4 先验信息的不同情景Fig. 4 Priori information of three scenarios

2.2 模型构建

本研究使用R包simmr来拟合所有的同位素贝叶斯混合模型(iter=50000,burn=1000,thin=10,n.chain=4),并使用JAGS(Just Another Gibbs Sampler)从后验分布中提取样本。通过马氏链(Markov chain,MCMC)轨迹图检验法、Geweke检验法和Gelman检验法等多种方法诊断了马氏链的收敛性,进一步估计有关参数或者进行其他统计推断。模型结果表明(表 1),蒙古鲌鱼食性功能群的浮游饵料鱼的比重最高,其次为底栖饵料鱼; 蒙古鲌浮游动物食性功能群的浮游动物和浮游饵料鱼的比重高; 蒙古鲌底栖动物食性功能群的底栖动物能量贡献比例最高,结果整体符合预期。

表1 模型结果Tab. 1 Model results

3 模型性能评价

3.1 整体性能

在模型选择阶段,常见指标为偏差信息量准则(Deviance information criterion,DIC)。DIC是等级模型化的赤池信息量准则(Akaike information criterion,AIC),被广泛应用于由马尔可夫链蒙特卡洛(MCMC)模拟出的后验分布的贝叶斯模型选择问题[36]。和赤池信息量准则一样,偏差信息量准则是随样本容量增加的渐近近似,只应用于后验分布呈多元正态分布的情况[11,37]。一般而言,偏差信息量准则的值越小,模型越好。这一准则的优点是它很容易从马尔可夫链蒙特卡洛(MCMC)模拟产生的样本中计算出来[36]。在一组供选择的模型中,如果观测值样本数不同的时候,可使用校准的DIC(DICcor)进行比较(公式3)[11,37],公式如下:

式中,N表示每个模型的样本数,n为比较的模型中样本数最低的数据。

例如,本研究中蒙古鲌鱼食性功能群、浮游动物食性功能群和底栖动物食性功能群的3个模型的样本数分别为11、10和7; 因此,计算DICcor的时候,取值n=7。由此模型整体性能评价结果表明(表 2),蒙古鲌鱼食性功能群、浮游动物食性功能群和底栖动物食性功能群的拟合总体较好(90% Coverage均为100%,90% Coverage为观测值落入90%后验分布中的比例)。但3个模型比较看,鱼食性功能群的拟合模型性能相对较差(DIC和DICcor最高)。需要特别指出的,在一组供选择的模型中,最优化的模型的选择都是具有相对性的,并不是说所选择的模型就一定足够精确。

表2 模型整体性能评价Tab. 2 Model performance

3.2 预测值与观测值

针对消费者营养来源研究,模型性能核心是指营养来源比例在不同方法下的预测效果,但是实践中很难获得消费者营养来源贡献的合理观测值。因此,通过评估消费者同位素的观测值(y)与预测值之间的匹配程度,可以间接评估不同食物来源比例的预测效果[41]。消费者稳定同位素预测值可通过模型中各种营养来源贡献比例和营养来源的稳定同位素来计算获得(公式4)[42]:

式中,n表示消费者使用的第n种食物来源,k是食物来源的数量(本研究中k=4),f代表通过食物来源分配方法预测的食物来源对消费者的比例贡献。δX表示食物来源中测得的同位素组成。

结合模型的特点,通常计算一个或多个基于预测值与测量值的指标来进行评价(表 3)。例如,均方根误差(RMSE),总结了测得值和预测值之间的平均差异,用来评估预测值的准确程度; 较低的RMSE值表明该模型具有较小的误差和更准确的预测。预测优度统计数据(G)也用于衡量模型的有效性。G值为1表示理想的预测。G值越接近1,模型的可靠性越高。G值为负表示该模型不太可靠。

表3 基于预测值与测量值的模型评价方法Tab. 3 Evaluation methods for the difference between model predicts and observations.

衡量模型的拟合程度(模型质量好坏),没有固定的标准。例如,MAE和RMSE一样,衡量的是真实值与预测值的偏离的绝对大小情况,需要结合真实值的量纲才能判断差异; 而MAPE衡量的是偏离的相对大小(即百分率)。相对来说,MAE和MAPE不容易受极端值的影响; 而MSE和RMSE采用误差的平方,会放大预测误差,所以对于离群数据更敏感。MAPE使用百分率来衡量偏离的大小,容易理解和解读。而MAE/RMSE需要结合真实值的量纲才能判断差异。

本研究的模型评价结果表明(表 4),蒙古鲌鱼食性功能群、浮游动物食性功能群和底栖动物食性功能群的拟合总体较好(G>0.8),但3个模型比较看,鱼食性功能群的拟合模型性能相对较差。不管用哪个指标,评估模型的好坏都不能够脱离具体的应用场景和具体的数据集。单纯地评判哪个模型好坏,是基本上没有意义的。

表4 模型预测值与观测值差异评价Tab. 4 Evaluation of the difference between model predicts and observations

3.3 先验信息与后验信息

贝叶斯混合模型的主要统计数据是消费者营养来源的概率分布,而不是消费者同位素值的预测,因此评估混合模型对消费者营养来源贡献估计的后验信息及其受先验信息的影响(图 5)。后验信息与先验信息的本质是营养贡献量的概率分布,因此可通过信息理论的相关概念和方法进行评价(表 5)。

图5 消费者不同营养功能群营养来源贡献的先验信息与后验信息Fig. 5 Priori information and posteriori information of trophic groups

表5 主要信息理论度量与统计Tab. 5 The important information theory measures and statistics

信息量用来度量一个信息的多少,而信息熵(Shannon Entropy)则用来描述一个来源信息的不确定度,也是信息来源的信息量期望[44]。基于信息理论的方法,可以衡量不同营养来源概率分布之间的先验信息与后验信息的异同。进一步可评估营养来源数据的改变、消费者功能群的改变和营养来源强度的改变等对消费者营养来源信息确定和解释的影响。例如,相对熵可以衡量先验信息与后验信息之间的距离,当先验信息与后验信息分布相同时,它们的相对熵为零[43]。在混合模型中,先验信息通常较后验信息丰富,因此当计算KL散度(Kullback-Leibler divergence),可获得比最大熵更高的信息增益,即后验与先验信息相差很大时,就会出现较高的值[43],而信息熵则代表了信息增益的有限上界。KL散度(信息增益)增高由后验信息与先验信息不同的均值和置信区间的差异变大决定(后验信息不支持先验信息或者后验信息强烈支持先验信息但置信区间较窄),因此需结合分布信息图一同解释(图 5)。

本研究依据主要信息理论度量与统计方法,评估了鲌贝叶斯混合模型的先验信息和后验信息结果(表 6),先验信息和后验信息不同的信息熵值、不为0的相对熵值及较大的交叉熵值表明鲌贝叶斯混合模型的先验信息和后验信息存在一定程度差异。首先,估计每个信息源的边际(Marginal)贡献所占的比例,以反映文献中混合模型是如何解释的(例如大多数混合模型程序的输出包括边际均值和可信度区间)。营养来源1%水平贡献信息测量应在0.01的离散区间内计算。其次,通过所有来源贡献的联合后验分布计算的,说明了关于来源相互间有条件依赖的信息增益,如食性组成权衡,及边际贡献的变化。我们通过比较食性组成贡献的先验分布和后验分布的样本来计算总体测量值,并使用The isometric logratio transformation等轴测对数比进行转换[48]。等轴测对数比率将比例转换为变量,这些变量根据多元正态分布近似分布,使我们能够使用2个多元正态分布散度的分析方程来计算信息增益。从概念上来说,我们可以把等轴测对数比变换看作是将食性组成比例延伸到一个可能有无限边界的坐标空间。等长对数比变换意味着KL散度可以无限增加,食性组成比例的变化非常微小。因此,我们在log2尺度上用KL散度来表示信息增益[11,37]。

表6 先验信息与后验信息的信息理论统计结果Tab. 6 The important information theory measures and statistics

此外,度量2个概率分布之间的距离,还有一些其他通用的散度指标可供参考(表 7)。例如Hellinger distance[49]可量化先验分布和后验分布之间差异。如果先验分布和后验分布在其参数空间中具有相同的密度,则为零。如果一个分布在任何地方都是零密度,而另一个分布是正密度,那么Hellinger distance为1[45]。

表7 概率分布相似度评价的几种距离法Tab. 7 Some distance methods to evaluate similarity of probabilistic distribution

4 结论

本文综述了在拟合和评估贝叶斯混合模型时遵循的最佳实践(图 6),及如何直接避免同位素构建消费者营养溯源分析中的诸多技术问题[9]。种群营养功能群的识别,有助于减小消费者同位素变异所带来的不确定性[19,20]; 营养来源矫正可以避免营养来源和消费者同位素之间的共线性,及不同营养来源之间的差异水平等情况[27]; 数据质量检验可以帮助剔除异常数据、提高数据总体质量、检验数据质量是否符合建模需要[28,29]。动态的生态系统与消费者相对快速的能量需求调整,可强烈地影响食物网结构与功能[50]。因此,在开展食物网研究前,需要投入大量的时间和精力进行数据收集和分析,了解野外调查过程中哪些数据满足了前期的假定条件,以确保建模分析的准确性[32]。在尝试数据建模之前,必须考虑环境中潜在食物来源的丰度、生物量、消费者的摄食行为习性、胃肠含物等重要的先验信息[33]。完整重现模型选择和模型评价的过程,是建模、训练、验证、评价的必要条件。模型选择是根据一组不同复杂度的模型表现,从中挑选最好的模型; 模型评价则是在选择模型后,评价其预测误差[11,37]。根据具体研究,在实践过程中,评价指标多种多样,且分别刻画了相对“真实模型”的信息损失。由于真实模型的未知性,这些评价仅反应现有模型构建过程中相对较好的性能,所以具体问题仍需具体分析。

图6 模型构建与性能评估流程图Fig. 6 Fishbone diagram of quantifying the quality of model construction and prediction

通过侧重于模型预测摄食者同位素值能力的评估方法,可以判断模型的拟合质量。此外,鉴于贝叶斯模型的特点,即如果先验信息误差较低,则信息混合模型的后验分布趋同于先验信息,进一步进行了基于信息理论和概率距离统计方法的评估,为同位素混合模型的输出结果的质量提供互补的评估方法[37]。这些方法的综合运用,将进一步提高消费者营养溯源准确性,为更深刻地认识食物网规律提供科学支撑[11,37]。