基于極大似然法的土壤重金屬刪失數(shù)據(jù)的相關(guān)性

馮 旭,孫大荃,李仁英,汪麗軍,黃利東* (.南京信息工程大學(xué)農(nóng)業(yè)資源與環(huán)境系,江蘇 南京 0044；.捷克科學(xué)院水土研究中心,捷克 布杰約維采 7005；.內(nèi)蒙古科爾沁右翼前旗農(nóng)牧業(yè)科學(xué)技術(shù)發(fā)展中心,內(nèi)蒙古 科右前旗 770)

土壤環(huán)境中的痕量重金屬作為土壤污染物或資源一直備受關(guān)注[1-9].由于痕量物質(zhì)的濃度較低,經(jīng)常發(fā)生樣品濃度低于檢測(cè)限的情況[10],此類情況下的樣本濃度數(shù)據(jù)被稱作左刪失數(shù)據(jù),在此情形下,研究者難以獲取數(shù)據(jù)的完整信息,給后續(xù)的數(shù)據(jù)分析工作帶來(lái)困難.針對(duì)一維刪失數(shù)據(jù),相關(guān)研究相繼利用參數(shù)模型[11-12]、半?yún)?shù)模型[13-14]以及非參數(shù)模型[15]對(duì)參數(shù)如均值和方差進(jìn)行了有效估計(jì)[16].然而研究測(cè)定的土壤痕量物質(zhì)常包含多個(gè)指標(biāo),且指標(biāo)間的相關(guān)性有著重要的研究?jī)r(jià)值.比如科研者不僅關(guān)注重金屬的濃度問(wèn)題,而且也研究重金屬的同源性問(wèn)題,以及修復(fù)措施對(duì)多種重金屬濃度協(xié)同變化的影響,這些信息都與變量之間的相關(guān)性密切相關(guān).如果二維數(shù)據(jù)都有刪失情況發(fā)生,那么相關(guān)性的估計(jì)將面臨更大挑戰(zhàn).在實(shí)際中,研究人員可能替換或刪除掉刪失的樣本以估計(jì)相關(guān)系數(shù),這樣做可能會(huì)產(chǎn)生有偏的結(jié)果.另外,剔除刪失樣本將會(huì)造成數(shù)據(jù)信息的浪費(fèi),達(dá)不到信息合理利用的目的.

痕量物質(zhì)的樣品濃度通常成對(duì)數(shù)正態(tài)分布[17-18],目前對(duì)符合對(duì)數(shù)正態(tài)分布的刪失數(shù)據(jù)之間相關(guān)性的報(bào)道較少.因此,本研究基于二維對(duì)數(shù)正態(tài)分布,考慮二維刪失的不同情況下的似然函數(shù),利用極大似然法(MLE)對(duì)相關(guān)系數(shù)進(jìn)行估計(jì),以澳大利亞土壤普查數(shù)據(jù)中的 Ag、Hg、Te、Hf為例(都有刪失)進(jìn)行方法示范和應(yīng)用.

1 材料與方法

1.1 分布形式與似然函數(shù)

二維對(duì)數(shù)正態(tài)分布的概率密度公式見(jiàn)式(1)[19]:

式(1)中呈二維對(duì)數(shù)正態(tài)分布的兩組數(shù)據(jù)分別記作x1與x2,y1=lnx1,y2=lnx2.μx1、μx2、σx1、σx2x2分別為x1、x2的均值和標(biāo)準(zhǔn)差,μy1、μy2、σy1、σy2分別為y1、y2的均值和標(biāo)準(zhǔn)差,μy1、μy2、σy1、σy2通過(guò)一維刪失對(duì)數(shù)正態(tài)分布的MLE進(jìn)行估計(jì),ρy1y2為y1、y2的相關(guān)系數(shù).

將二維對(duì)數(shù)正態(tài)數(shù)據(jù)的刪失情況分為4種:第1種情況x1、x2皆未刪失;第2種情況x1刪失x2未刪失;第3種情況x2刪失x1未刪失;第4種情況x1、x2皆刪失.似然函數(shù) L(θ)如下:

式(2)中n為樣本容量.Lx1模擬對(duì)x1檢測(cè)方法的檢測(cè)限,Lx2模擬對(duì)x2檢測(cè)方法的檢測(cè)限.

1.2 模擬研究的方法

探究樣本容量、相關(guān)系數(shù)、刪失比例和干擾項(xiàng)等4個(gè)因素對(duì)MLE估計(jì)值準(zhǔn)確性的影響,與刪除法、替換法對(duì)比,檢驗(yàn)MLE的估計(jì)效果.

1.2.1 樣本容量對(duì)MLE相關(guān)系數(shù)估計(jì)值準(zhǔn)確性的影響 通過(guò)生成不同數(shù)量的二維對(duì)數(shù)正態(tài)隨機(jī)數(shù),研究樣本容量大小對(duì)相關(guān)系數(shù)估計(jì)的準(zhǔn)確性的影響,確定本研究最佳樣本容量.本研究對(duì)不同均值參數(shù)生成的隨機(jī)數(shù)進(jìn)行了模擬研究,估計(jì)結(jié)果基本一致,因此文中以一組均值為3的隨機(jī)數(shù)為例,使,避免模擬計(jì)算中出現(xiàn)大量含小數(shù)或大數(shù)值運(yùn)算,提高運(yùn)算效率,使,保證隨機(jī)數(shù)有較大的變化幅度,方便對(duì)本文方法的客觀評(píng)價(jià).上述參數(shù)選擇相同均值和方差,有利于區(qū)分估計(jì)結(jié)果準(zhǔn)確性變化的主導(dǎo)來(lái)源(一維參數(shù)估計(jì)的準(zhǔn)確性或似然函數(shù)).將每個(gè)樣本容量下的不同刪失比例的數(shù)據(jù)的估計(jì)值整合在一起進(jìn)行對(duì)比分析.

在上述的基礎(chǔ)上固定樣本容量,通過(guò)設(shè)置Lx1與Lx2的大小調(diào)控刪失比例,并設(shè)置不同參數(shù)隨機(jī)數(shù)下的模擬試驗(yàn)并進(jìn)行比較,排除相關(guān)系數(shù)估計(jì)值受樣本大小等因素的干擾,詳細(xì)研究刪失比例對(duì)相關(guān)系數(shù)估計(jì)值(ry1y2)準(zhǔn)確性的影響,提供更全面客觀的模擬結(jié)果.

1.2.2 總體相關(guān)系數(shù)對(duì)MLE相關(guān)系數(shù)估計(jì)值準(zhǔn)確性的影響 在上述的基礎(chǔ)上,通過(guò)控制協(xié)方差,調(diào)控總體相關(guān)系數(shù)的大小,重復(fù)模擬 1000次,給出隨刪失比例增長(zhǎng)估計(jì)值的變化區(qū)間,由于估計(jì)結(jié)果具有對(duì)稱性,固定x2的刪失比例為50%,控制x1的刪失比例研究總體相關(guān)系數(shù)變化(-1～1,間隔為0.1)對(duì)估計(jì)值準(zhǔn)確性的影響.隨機(jī)數(shù)的相關(guān)系數(shù)為-1～1均勻分布,因此以總體相關(guān)系數(shù)為0.5、-0.5為例,對(duì)比不同總體相關(guān)系數(shù)下,估計(jì)值準(zhǔn)確性隨刪失比例變化的差異.

1.2.3 干擾項(xiàng)對(duì)MLE相關(guān)系數(shù)估計(jì)值準(zhǔn)確性的影響 在上述的基礎(chǔ)上,通過(guò)樣本數(shù)據(jù)加正態(tài)隨機(jī)數(shù)(均值為0,方差為10%～50%LOD,間隔為20%LOD),引入干擾項(xiàng),研究干擾項(xiàng)對(duì)MLE估計(jì)值準(zhǔn)確性的影響,測(cè)試本文方法的魯棒性.

1.2.4 不同方法估計(jì)刪失數(shù)據(jù)相關(guān)系數(shù)的準(zhǔn)確性比較 實(shí)際樣品測(cè)定常會(huì)受到方法或儀器的影響而使測(cè)定結(jié)果有一定的變動(dòng)范圍,因此通過(guò)多次取樣(1000次),比較替換法、刪除法和本文方法的準(zhǔn)確性和精確性.

1.3 評(píng)價(jià)方法

準(zhǔn)確性:合并同樣本容量下各刪失比例數(shù)據(jù)的相關(guān)系數(shù)估計(jì)值并求其均值,從而比較不同樣本容量下相關(guān)系數(shù)估計(jì)值準(zhǔn)確性的整體水平;通過(guò)估計(jì)值的等值線分布與估計(jì)值取值范圍,評(píng)價(jià)與探究估計(jì)值的準(zhǔn)確性與偏差的趨勢(shì).

精確性:以琴型圖的分位數(shù)間隔與散點(diǎn)分布情況表征估計(jì)結(jié)果的精確性.

魯棒性:通過(guò)bootstrap重采樣,對(duì)未經(jīng)刪失的數(shù)據(jù)進(jìn)行 95%置信水平的區(qū)間估計(jì),評(píng)價(jià)估計(jì)結(jié)果的可信程度,分析加干擾項(xiàng)后本文MLE是否出現(xiàn)更多超出置信區(qū)間的估計(jì)值,檢驗(yàn)方法的魯棒性.

1.4 實(shí)例分析

利用澳大利亞土壤普查[20]中的 Ag、Hg、Te、Hf 等(http://dx.doi.org/10.11636/Record.2011.020)痕量物質(zhì)(呈對(duì)數(shù)正態(tài)分布[17-18],且都存在不同程度的刪失)對(duì)本方法進(jìn)行評(píng)價(jià).估算 4種元素之間的相關(guān)性,對(duì)比刪除法、替換法與MLE的估計(jì)效果.

1.5 分析方法

數(shù)據(jù)分析通過(guò) R語(yǔ)言(4.1.2版本)實(shí)現(xiàn),隨機(jī)數(shù)由 compositions(2.0-2)中的 rlnorm.rplus()生成,最大似然估計(jì)通過(guò)maxLik (1.5-2)實(shí)現(xiàn),替換法和刪除法通過(guò) stats (3.6.2)中的 cor()計(jì)算,繪圖工具選用OriginPro 2021與R語(yǔ)言中的ggplot2 (3.3.5).

2 結(jié)果與討論

2.1 樣本容量對(duì)MLE相關(guān)系數(shù)估計(jì)值準(zhǔn)確性的影響

由圖1可見(jiàn),n=50時(shí),各刪失比例(10%～90%,間隔為10%)共100個(gè)相關(guān)系數(shù)估計(jì)值,如圖中均值線所示,均值()在0.37左右,隨著n的增大,相關(guān)系數(shù)估計(jì)值逐漸聚集,逐漸向設(shè)定的 ρy1y2靠近,當(dāng)n=2000時(shí)均值變化較小,估計(jì)效果趨于穩(wěn)定,當(dāng)n進(jìn)一步增大時(shí),估計(jì)結(jié)果進(jìn)一步收斂.

圖1 樣本容量對(duì)MLE估計(jì)準(zhǔn)確性的影響Fig.1 The influence of sample size on the accuracy of MLE estimates

從樣本攜帶信息的角度,對(duì)于樣本數(shù)量較小的數(shù)據(jù),攜帶的信息量較少,隨著刪失比例的增長(zhǎng)使Fisher信息[21-22]進(jìn)一步減少,對(duì)估計(jì)結(jié)果的均值產(chǎn)生較大影響,因此當(dāng)樣本容量為50時(shí),本方法的估計(jì)效果較差,但經(jīng)驗(yàn)證本文提出的方法,盡管在樣本容量為50時(shí)仍優(yōu)于刪失法與替換法.隨著樣本數(shù)量的增多,偏差逐漸減小[23-24].對(duì)于較大樣本數(shù)據(jù)的估計(jì),MLE達(dá)到了預(yù)期的效果[25-26].整體來(lái)看,本研究的 MLE具有漸進(jìn)無(wú)偏性與一致性.樣本容量達(dá)到2000時(shí),估計(jì)結(jié)果的均值與方差已基本收斂,故后續(xù)模擬研究的樣本量設(shè)為2000.

2.2 刪失比例對(duì)MLE相關(guān)系數(shù)估計(jì)值準(zhǔn)確性的影響

為保證模擬實(shí)驗(yàn)的客觀性,以不同的隨機(jī)數(shù)組合(表1)的3次模擬實(shí)驗(yàn)為例(圖2)研究刪失比例對(duì)似然估計(jì)值準(zhǔn)確性的影響.不同隨機(jī)樣本下相關(guān)系數(shù)的似然估計(jì)值隨刪失比例變化情況基本一致.刪失比例在 0～35%左右時(shí),估計(jì)值相對(duì)于真實(shí)值偏高,當(dāng)刪失比例到達(dá) 60%左右,等高線密度增加,相關(guān)系數(shù)的估計(jì)值變化加快.隨刪失比例的上升,相關(guān)系數(shù)的估計(jì)值逐漸減小,估計(jì)結(jié)果的分布趨勢(shì)大致呈不規(guī)則的對(duì)稱,在 ρy1y2取不同值時(shí),有著相似的變化趨勢(shì),但隨 ρy1y2變大,刪失比例高于 80%的估計(jì)值有較大的偏差出現(xiàn).

表1 干擾項(xiàng)對(duì)MLE方法的魯棒性檢驗(yàn)Table 1 Robustness test of MLE method by introducing disturbance term

隨著刪失比例的上升,似然估計(jì)值的變化較為均勻,且不同隨機(jī)數(shù)樣本(總體相關(guān)系數(shù)一致)的相關(guān)系數(shù)估計(jì)結(jié)果相似(圖2),只有標(biāo)準(zhǔn)差達(dá)到3時(shí),等高線的線條相對(duì)趨于平緩,差異性較小,說(shuō)明隨機(jī)數(shù)在滿足分布的前提下,本身的大小與離散程度并不會(huì)對(duì)估計(jì)結(jié)果的準(zhǔn)確性造成較大影響.而與預(yù)期中不同的是,估計(jì)值的等高線圖并不是嚴(yán)格對(duì)稱的,這是由于與常規(guī)的相關(guān)系數(shù)計(jì)算不同,極大似然法基于似然函數(shù)得出最優(yōu)解,因此估計(jì)值可能會(huì)有輕微變化與波動(dòng),但偏差控制在0.06以內(nèi),并未對(duì)參數(shù)的估計(jì)造成較大的影響[27].

圖2 刪失比例對(duì)MLE相關(guān)系數(shù)估計(jì)值準(zhǔn)確性的影響Fig.2 Influence of percent censored on the accuracy of MLE correlation coefficient estimation

2.3 總體相關(guān)系數(shù)對(duì)MLE估計(jì)值準(zhǔn)確性的影響

如圖3所示,總體相關(guān)系數(shù)取0.5時(shí),似然函數(shù)的集散程度受x2的影響較大,在刪失比例較低處,1000次重復(fù)的相關(guān)系數(shù)的估計(jì)值較為集中,當(dāng)刪失比例到達(dá)60%左右,估計(jì)值逐漸趨于離散.總體相關(guān)系數(shù)取-0.5時(shí),估計(jì)值的集散有著相同的趨勢(shì),而對(duì)比總體相關(guān)系數(shù)為0.5的估計(jì)值,分布更為集中,離散程度變化更加明顯.

圖3 不同刪失比例下對(duì)MLE估計(jì)值準(zhǔn)確性與總體相關(guān)系數(shù)的影響Fig.3 Change of accuracy of MLE estimates with population correlation under different censored percentage

可以看出相關(guān)系數(shù)的估計(jì)值受總體相關(guān)水平的影響,類似的研究中也提到了這一現(xiàn)象[28],可能是隨著刪失比例的變化,對(duì)樣本的均值和方差的估計(jì)造成了影響,由皮爾遜系數(shù)的計(jì)算公式可知,總體相關(guān)系數(shù)與協(xié)方差的大小密切相關(guān),總體相關(guān)系數(shù)的變化使依據(jù)協(xié)方差所生成隨機(jī)數(shù)的統(tǒng)計(jì)量出現(xiàn)一定差異,導(dǎo)致估計(jì)方法在不同的總體相關(guān)水平下對(duì)刪失數(shù)據(jù)相關(guān)系數(shù)的估計(jì)效果出現(xiàn)變化.基于此猜想,對(duì)均值與標(biāo)準(zhǔn)差的估計(jì)值進(jìn)行了檢驗(yàn),總體相關(guān)系數(shù)的確對(duì)其估計(jì)效果造成了影響,在均值與方差的估計(jì)值偏差較大處,相關(guān)系數(shù)的估計(jì)值也隨之出現(xiàn)了較大的變化,可見(jiàn)方法對(duì)統(tǒng)計(jì)量估計(jì)值的準(zhǔn)確性有一定的依賴性,因此對(duì)于一維統(tǒng)計(jì)量的估計(jì)方法仍需改進(jìn).但整體而言,MLE對(duì)均值與方差估計(jì)較為準(zhǔn)確[29],相關(guān)系數(shù)的估計(jì)值受到的影響不大.

2.4 干擾項(xiàng)對(duì)MLE的魯棒性檢驗(yàn)

實(shí)際應(yīng)用中,檢測(cè)儀器和實(shí)驗(yàn)方法的誤差,可能會(huì)使得到的樣品數(shù)據(jù)出現(xiàn)偏差,添加干擾項(xiàng)對(duì)此類情況進(jìn)行了模擬.如上述結(jié)果表明,正態(tài)干擾項(xiàng)的施加并未對(duì)本研究的MLE造成較大影響,可看出方法具有較好的魯棒性[30].

2.5 不同方法相關(guān)系數(shù)估計(jì)值的精準(zhǔn)性比較

模擬過(guò)程重復(fù) 1000次后,相關(guān)系數(shù)如圖4所示[28].LOD替換法與LOD/2替換法有著相似的變化趨勢(shì),估計(jì)值隨刪失比例的提高逐漸降低,當(dāng)下降到0.3左右,變化趨于穩(wěn)定.刪除法在2組數(shù)據(jù)的刪失比例較低時(shí)便無(wú)法進(jìn)行有效估計(jì),相關(guān)系數(shù)估計(jì)值隨刪失比例的增加而迅速減小并逐漸分散.MLE方法的估計(jì)結(jié)果變化較小,估計(jì)值圍繞在總體相關(guān)系數(shù)的周圍,基本聚集在0.4～0.6,估計(jì)效果最好.當(dāng)刪除法與替換法的 ρy1y2絕對(duì)值較大時(shí),變化尤其明顯,而MLE受到的影響較小,可以保證在 ρy1y2變化的情況下,數(shù)據(jù)刪失比例達(dá)到80%,仍可以進(jìn)行有效估計(jì).

圖4 不同方法估計(jì)刪失數(shù)據(jù)相關(guān)系數(shù)的準(zhǔn)確性比較Fig.4 Comparison of accuracy of correlation coefficient estimated by different methods

經(jīng)重復(fù),替換法估計(jì)值不斷下降后趨于穩(wěn)定,可能是因?yàn)樘鎿Q導(dǎo)致的樣本整體均值水平升高,使估計(jì)值不斷降低,當(dāng)刪失比例達(dá)到一定程度,替換的固定值占樣本數(shù)據(jù)的絕大部分,相關(guān)系數(shù)趨于穩(wěn)定.刪除法由于刪失部分樣本信息的直接丟失,且未有替換法中的固定值填補(bǔ),導(dǎo)致估計(jì)值的不斷減小.替換法優(yōu)于刪除法.本研究的 MLE估計(jì)值較為準(zhǔn)確且穩(wěn)定,優(yōu)于替換法和刪除法,有較高的精準(zhǔn)度[31].

2.6 實(shí)例分析

以澳大利亞土壤普查數(shù)據(jù)中的Ag,Hg,Te,Hf(刪失比例分別約為20%、30%、40%、50%)為例,經(jīng)檢驗(yàn)樣本數(shù)據(jù)皆符合對(duì)數(shù)正態(tài)分布,刪失情況如表2.

表2 2011年澳大利亞土壤普查數(shù)據(jù)的樣本容量與具體刪失情況Table 2 Sample size and censored percentage of soil survey data in Australia in 2011

如圖5所示,除去刪除法,替換法和MLE對(duì)Hg與Te,Hg與Hf以及Ag和Hg相關(guān)系數(shù)的估計(jì)值較為集中,而Hf與Ag以及Hf和Te的相關(guān)系數(shù)的估計(jì)值較為分散,MLE的估計(jì)值常大于刪除法與替換法.

圖5 不同方法對(duì)土壤重金屬(含刪失)相關(guān)性估計(jì)比較Fig.5 Comparison of correlation coefficient of soil heavy metals (including censored) estimated by different methods

4種元素的樣本容量達(dá)到了5226(表2),各種方法的估計(jì)結(jié)果已大致收斂,因此除了刪除法,替換法與 MLE的結(jié)果有著較為相似的趨勢(shì)(圖5)[32].刪除法明顯受到了刪失比例的影響,在對(duì)含有刪失比例達(dá)到51.5%的Hf元素的相關(guān)系數(shù)估計(jì)中,刪除法與其他方法的估計(jì)值相差較大,而 LOD替換法與LOD/2替換法估計(jì)值相對(duì)于MLE偏低[33],尤其是在對(duì)相關(guān)系數(shù)較高與刪失比例較高的元素進(jìn)行估計(jì)時(shí)較為明顯,這與模擬研究中的結(jié)果基本吻合.通過(guò)相關(guān)系數(shù)估計(jì),Ag與Hg有著較高的相關(guān)性,或許有相同的來(lái)源[34],Hg與Hf相關(guān)系數(shù)趨近于0,基本沒(méi)有關(guān)聯(lián).

本文算法和代碼是開(kāi)放的提出的方法僅對(duì)左刪失數(shù)據(jù)進(jìn)行了針對(duì)性研究,但在理論上廣泛適用于各種刪失機(jī)制,可以實(shí)現(xiàn)二維刪失數(shù)據(jù)的相關(guān)性無(wú)偏估計(jì),為了解土壤痕量物質(zhì)(不局限于重金屬)的關(guān)聯(lián)性提供基礎(chǔ).

3 結(jié)論

3.1 樣本容量越大,MLE的結(jié)果越準(zhǔn)確,當(dāng)達(dá)到一定樣本容量(2000)時(shí)估計(jì)效果趨于穩(wěn)定.

3.2 隨著刪失比例增加,基于對(duì)數(shù)正態(tài)分布的MLE的相關(guān)系數(shù)變化較小,且總體相關(guān)系數(shù)的變化未對(duì)其造成較大影響,表明本方法的穩(wěn)定性和一致性.

3.3 隨機(jī)干擾項(xiàng)的引入未對(duì)估計(jì)結(jié)果造成較大影響,表明本文提出的MLE有較強(qiáng)的魯棒性.

3.4 隨著刪失比例的提升,刪除法和替換法的結(jié)果準(zhǔn)確性變差,而MLE準(zhǔn)確性明顯優(yōu)于上述兩種方法.

3.5 實(shí)際數(shù)據(jù)驗(yàn)證表明,本文方法相比于替換法和刪除法,擬合結(jié)果更加穩(wěn)定,受刪失比例變化的影響較小,Ag與Hg具有較高的相關(guān)性,可能具有相同來(lái)源,Hg與Hf之間相關(guān)系數(shù)趨近于0,基本無(wú)關(guān)聯(lián).