基于极大似然法的土壤重金属删失数据的相关性

冯 旭,孙大荃,李仁英,汪丽军,黄利东* (.南京信息工程大学农业资源与环境系,江苏 南京 0044；
.捷克科学院水土研究中心,捷克 布杰约维采 7005；
.内蒙古科尔沁右翼前旗农牧业科学技术发展中心,内蒙古 科右前旗 770)

土壤环境中的痕量重金属作为土壤污染物或资源一直备受关注[1-9].由于痕量物质的浓度较低,经常发生样品浓度低于检测限的情况[10],此类情况下的样本浓度数据被称作左删失数据,在此情形下,研究者难以获取数据的完整信息,给后续的数据分析工作带来困难.针对一维删失数据,相关研究相继利用参数模型[11-12]、半参数模型[13-14]以及非参数模型[15]对参数如均值和方差进行了有效估计[16].然而研究测定的土壤痕量物质常包含多个指标,且指标间的相关性有着重要的研究价值.比如科研者不仅关注重金属的浓度问题,而且也研究重金属的同源性问题,以及修复措施对多种重金属浓度协同变化的影响,这些信息都与变量之间的相关性密切相关.如果二维数据都有删失情况发生,那么相关性的估计将面临更大挑战.在实际中,研究人员可能替换或删除掉删失的样本以估计相关系数,这样做可能会产生有偏的结果.另外,剔除删失样本将会造成数据信息的浪费,达不到信息合理利用的目的.

痕量物质的样品浓度通常成对数正态分布[17-18],目前对符合对数正态分布的删失数据之间相关性的报道较少.因此,本研究基于二维对数正态分布,考虑二维删失的不同情况下的似然函数,利用极大似然法(MLE)对相关系数进行估计,以澳大利亚土壤普查数据中的 Ag、Hg、Te、Hf为例(都有删失)进行方法示范和应用.

1.1 分布形式与似然函数

二维对数正态分布的概率密度公式见式(1)[19]:

式(1)中呈二维对数正态分布的两组数据分别记作x1与x2,y1=lnx1,y2=lnx2.μx1、μx2、σx1、σx2x2分别为x1、x2的均值和标准差,μy1、μy2、σy1、σy2分别为y1、y2的均值和标准差,μy1、μy2、σy1、σy2通过一维删失对数正态分布的MLE进行估计,ρy1y2为y1、y2的相关系数.

将二维对数正态数据的删失情况分为4种:第1种情况x1、x2皆未删失;第2种情况x1删失x2未删失;第3种情况x2删失x1未删失;第4种情况x1、x2皆删失.似然函数 L(θ)如下:

式(2)中n为样本容量.Lx1模拟对x1检测方法的检测限,Lx2模拟对x2检测方法的检测限.

1.2 模拟研究的方法

探究样本容量、相关系数、删失比例和干扰项等4个因素对MLE估计值准确性的影响,与删除法、替换法对比,检验MLE的估计效果.

1.2.1 样本容量对MLE相关系数估计值准确性的影响 通过生成不同数量的二维对数正态随机数,研究样本容量大小对相关系数估计的准确性的影响,确定本研究最佳样本容量.本研究对不同均值参数生成的随机数进行了模拟研究,估计结果基本一致,因此文中以一组均值为3的随机数为例,使,避免模拟计算中出现大量含小数或大数值运算,提高运算效率,使,保证随机数有较大的变化幅度,方便对本文方法的客观评价.上述参数选择相同均值和方差,有利于区分估计结果准确性变化的主导来源(一维参数估计的准确性或似然函数).将每个样本容量下的不同删失比例的数据的估计值整合在一起进行对比分析.

在上述的基础上固定样本容量,通过设置Lx1与Lx2的大小调控删失比例,并设置不同参数随机数下的模拟试验并进行比较,排除相关系数估计值受样本大小等因素的干扰,详细研究删失比例对相关系数估计值(ry1y2)准确性的影响,提供更全面客观的模拟结果.

1.2.2 总体相关系数对MLE相关系数估计值准确性的影响 在上述的基础上,通过控制协方差,调控总体相关系数的大小,重复模拟 1000次,给出随删失比例增长估计值的变化区间,由于估计结果具有对称性,固定x2的删失比例为50%,控制x1的删失比例研究总体相关系数变化(-1～1,间隔为0.1)对估计值准确性的影响.随机数的相关系数为-1～1均匀分布,因此以总体相关系数为0.5、-0.5为例,对比不同总体相关系数下,估计值准确性随删失比例变化的差异.

1.2.3 干扰项对MLE相关系数估计值准确性的影响 在上述的基础上,通过样本数据加正态随机数(均值为0,方差为10%～50%LOD,间隔为20%LOD),引入干扰项,研究干扰项对MLE估计值准确性的影响,测试本文方法的鲁棒性.

1.2.4 不同方法估计删失数据相关系数的准确性比较 实际样品测定常会受到方法或仪器的影响而使测定结果有一定的变动范围,因此通过多次取样(1000次),比较替换法、删除法和本文方法的准确性和精确性.

1.3 评价方法

准确性:合并同样本容量下各删失比例数据的相关系数估计值并求其均值,从而比较不同样本容量下相关系数估计值准确性的整体水平;通过估计值的等值线分布与估计值取值范围,评价与探究估计值的准确性与偏差的趋势.

精确性:以琴型图的分位数间隔与散点分布情况表征估计结果的精确性.

鲁棒性:通过bootstrap重采样,对未经删失的数据进行 95%置信水平的区间估计,评价估计结果的可信程度,分析加干扰项后本文MLE是否出现更多超出置信区间的估计值,检验方法的鲁棒性.

1.4 实例分析

利用澳大利亚土壤普查[20]中的 Ag、Hg、Te、Hf 等(http://dx.doi.org/10.11636/Record.2011.020)痕量物质(呈对数正态分布[17-18],且都存在不同程度的删失)对本方法进行评价.估算 4种元素之间的相关性,对比删除法、替换法与MLE的估计效果.

1.5 分析方法

数据分析通过 R语言(4.1.2版本)实现,随机数由 compositions(2.0-2)中的 rlnorm.rplus()生成,最大似然估计通过maxLik (1.5-2)实现,替换法和删除法通过 stats (3.6.2)中的 cor()计算,绘图工具选用OriginPro 2021与R语言中的ggplot2 (3.3.5).

2.1 样本容量对MLE相关系数估计值准确性的影响

由图1可见,n=50时,各删失比例(10%～90%,间隔为10%)共100个相关系数估计值,如图中均值线所示,均值()在0.37左右,随着n的增大,相关系数估计值逐渐聚集,逐渐向设定的 ρy1y2靠近,当n=2000时均值变化较小,估计效果趋于稳定,当n进一步增大时,估计结果进一步收敛.

图1 样本容量对MLE估计准确性的影响Fig.1 The influence of sample size on the accuracy of MLE estimates

从样本携带信息的角度,对于样本数量较小的数据,携带的信息量较少,随着删失比例的增长使Fisher信息[21-22]进一步减少,对估计结果的均值产生较大影响,因此当样本容量为50时,本方法的估计效果较差,但经验证本文提出的方法,尽管在样本容量为50时仍优于删失法与替换法.随着样本数量的增多,偏差逐渐减小[23-24].对于较大样本数据的估计,MLE达到了预期的效果[25-26].整体来看,本研究的 MLE具有渐进无偏性与一致性.样本容量达到2000时,估计结果的均值与方差已基本收敛,故后续模拟研究的样本量设为2000.

2.2 删失比例对MLE相关系数估计值准确性的影响

为保证模拟实验的客观性,以不同的随机数组合(表1)的3次模拟实验为例(图2)研究删失比例对似然估计值准确性的影响.不同随机样本下相关系数的似然估计值随删失比例变化情况基本一致.删失比例在 0～35%左右时,估计值相对于真实值偏高,当删失比例到达 60%左右,等高线密度增加,相关系数的估计值变化加快.随删失比例的上升,相关系数的估计值逐渐减小,估计结果的分布趋势大致呈不规则的对称,在 ρy1y2取不同值时,有着相似的变化趋势,但随 ρy1y2变大,删失比例高于 80%的估计值有较大的偏差出现.

表1 干扰项对MLE方法的鲁棒性检验Table 1 Robustness test of MLE method by introducing disturbance term

随着删失比例的上升,似然估计值的变化较为均匀,且不同随机数样本(总体相关系数一致)的相关系数估计结果相似(图2),只有标准差达到3时,等高线的线条相对趋于平缓,差异性较小,说明随机数在满足分布的前提下,本身的大小与离散程度并不会对估计结果的准确性造成较大影响.而与预期中不同的是,估计值的等高线图并不是严格对称的,这是由于与常规的相关系数计算不同,极大似然法基于似然函数得出最优解,因此估计值可能会有轻微变化与波动,但偏差控制在0.06以内,并未对参数的估计造成较大的影响[27].

图2 删失比例对MLE相关系数估计值准确性的影响Fig.2 Influence of percent censored on the accuracy of MLE correlation coefficient estimation

2.3 总体相关系数对MLE估计值准确性的影响

如图3所示,总体相关系数取0.5时,似然函数的集散程度受x2的影响较大,在删失比例较低处,1000次重复的相关系数的估计值较为集中,当删失比例到达60%左右,估计值逐渐趋于离散.总体相关系数取-0.5时,估计值的集散有着相同的趋势,而对比总体相关系数为0.5的估计值,分布更为集中,离散程度变化更加明显.

图3 不同删失比例下对MLE估计值准确性与总体相关系数的影响Fig.3 Change of accuracy of MLE estimates with population correlation under different censored percentage

可以看出相关系数的估计值受总体相关水平的影响,类似的研究中也提到了这一现象[28],可能是随着删失比例的变化,对样本的均值和方差的估计造成了影响,由皮尔逊系数的计算公式可知,总体相关系数与协方差的大小密切相关,总体相关系数的变化使依据协方差所生成随机数的统计量出现一定差异,导致估计方法在不同的总体相关水平下对删失数据相关系数的估计效果出现变化.基于此猜想,对均值与标准差的估计值进行了检验,总体相关系数的确对其估计效果造成了影响,在均值与方差的估计值偏差较大处,相关系数的估计值也随之出现了较大的变化,可见方法对统计量估计值的准确性有一定的依赖性,因此对于一维统计量的估计方法仍需改进.但整体而言,MLE对均值与方差估计较为准确[29],相关系数的估计值受到的影响不大.

2.4 干扰项对MLE的鲁棒性检验

实际应用中,检测仪器和实验方法的误差,可能会使得到的样品数据出现偏差,添加干扰项对此类情况进行了模拟.如上述结果表明,正态干扰项的施加并未对本研究的MLE造成较大影响,可看出方法具有较好的鲁棒性[30].

2.5 不同方法相关系数估计值的精准性比较

模拟过程重复 1000次后,相关系数如图4所示[28].LOD替换法与LOD/2替换法有着相似的变化趋势,估计值随删失比例的提高逐渐降低,当下降到0.3左右,变化趋于稳定.删除法在2组数据的删失比例较低时便无法进行有效估计,相关系数估计值随删失比例的增加而迅速减小并逐渐分散.MLE方法的估计结果变化较小,估计值围绕在总体相关系数的周围,基本聚集在0.4～0.6,估计效果最好.当删除法与替换法的 ρy1y2绝对值较大时,变化尤其明显,而MLE受到的影响较小,可以保证在 ρy1y2变化的情况下,数据删失比例达到80%,仍可以进行有效估计.

图4 不同方法估计删失数据相关系数的准确性比较Fig.4 Comparison of accuracy of correlation coefficient estimated by different methods

经重复,替换法估计值不断下降后趋于稳定,可能是因为替换导致的样本整体均值水平升高,使估计值不断降低,当删失比例达到一定程度,替换的固定值占样本数据的绝大部分,相关系数趋于稳定.删除法由于删失部分样本信息的直接丢失,且未有替换法中的固定值填补,导致估计值的不断减小.替换法优于删除法.本研究的 MLE估计值较为准确且稳定,优于替换法和删除法,有较高的精准度[31].

2.6 实例分析

以澳大利亚土壤普查数据中的Ag,Hg,Te,Hf(删失比例分别约为20%、30%、40%、50%)为例,经检验样本数据皆符合对数正态分布,删失情况如表2.

表2 2011年澳大利亚土壤普查数据的样本容量与具体删失情况Table 2 Sample size and censored percentage of soil survey data in Australia in 2011

如图5所示,除去删除法,替换法和MLE对Hg与Te,Hg与Hf以及Ag和Hg相关系数的估计值较为集中,而Hf与Ag以及Hf和Te的相关系数的估计值较为分散,MLE的估计值常大于删除法与替换法.

图5 不同方法对土壤重金属(含删失)相关性估计比较Fig.5 Comparison of correlation coefficient of soil heavy metals (including censored) estimated by different methods

4种元素的样本容量达到了5226(表2),各种方法的估计结果已大致收敛,因此除了删除法,替换法与 MLE的结果有着较为相似的趋势(图5)[32].删除法明显受到了删失比例的影响,在对含有删失比例达到51.5%的Hf元素的相关系数估计中,删除法与其他方法的估计值相差较大,而 LOD替换法与LOD/2替换法估计值相对于MLE偏低[33],尤其是在对相关系数较高与删失比例较高的元素进行估计时较为明显,这与模拟研究中的结果基本吻合.通过相关系数估计,Ag与Hg有着较高的相关性,或许有相同的来源[34],Hg与Hf相关系数趋近于0,基本没有关联.

本文算法和代码是开放的提出的方法仅对左删失数据进行了针对性研究,但在理论上广泛适用于各种删失机制,可以实现二维删失数据的相关性无偏估计,为了解土壤痕量物质(不局限于重金属)的关联性提供基础.

3.1 样本容量越大,MLE的结果越准确,当达到一定样本容量(2000)时估计效果趋于稳定.

3.2 随着删失比例增加,基于对数正态分布的MLE的相关系数变化较小,且总体相关系数的变化未对其造成较大影响,表明本方法的稳定性和一致性.

3.3 随机干扰项的引入未对估计结果造成较大影响,表明本文提出的MLE有较强的鲁棒性.

3.4 随着删失比例的提升,删除法和替换法的结果准确性变差,而MLE准确性明显优于上述两种方法.

3.5 实际数据验证表明,本文方法相比于替换法和删除法,拟合结果更加稳定,受删失比例变化的影响较小,Ag与Hg具有较高的相关性,可能具有相同来源,Hg与Hf之间相关系数趋近于0,基本无关联.