基于極大似然法的土壤重金屬刪失數據的相關性

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

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

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

1 材料與方法

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 結果與討論

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 結論

3.1 樣本容量越大,MLE的結果越準確,當達到一定樣本容量(2000)時估計效果趨于穩定.

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

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

3.4 隨著刪失比例的提升,刪除法和替換法的結果準確性變差,而MLE準確性明顯優于上述兩種方法.

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