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

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

1 材料與方法

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

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

式(1)中呈二維對數(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通過一維刪失對數(shù)正態(tài)分布的MLE進(jìn)行估計,ρy1y2為y1、y2的相關(guān)系數(shù).

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

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

1.2 模擬研究的方法

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

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

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

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

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

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

1.3 評價方法

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

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

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

1.4 實例分析

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

1.5 分析方法

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

2 結(jié)果與討論

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

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

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

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

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

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

表1 干擾項對MLE方法的魯棒性檢驗Table 1 Robustness test of MLE method by introducing disturbance term

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

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

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

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

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

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

2.4 干擾項對MLE的魯棒性檢驗

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

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

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

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

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

2.6 實例分析

以澳大利亞土壤普查數(shù)據(jù)中的Ag,Hg,Te,Hf(刪失比例分別約為20%、30%、40%、50%)為例,經(jīng)檢驗樣本數(shù)據(jù)皆符合對數(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對Hg與Te,Hg與Hf以及Ag和Hg相關(guān)系數(shù)的估計值較為集中,而Hf與Ag以及Hf和Te的相關(guān)系數(shù)的估計值較為分散,MLE的估計值常大于刪除法與替換法.

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

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

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

3 結(jié)論

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

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

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

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

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