劉玥，施三支

（長春理工大學(xué) 理學(xué)院，長春 130022）

數(shù)據(jù)的缺失幾乎存在于所有的實證研究領(lǐng)域?？赡軙霈F(xiàn)數(shù)據(jù)缺失的原因有很多，包括調(diào)查無響應(yīng)、無法測量和單純的數(shù)據(jù)丟失等。但同時，現(xiàn)在各類研究不斷増多，這就需要不斷采集挖掘數(shù)據(jù)；數(shù)據(jù)量越大，數(shù)據(jù)缺失的可能性也就越大，缺失數(shù)據(jù)就會對統(tǒng)計推斷結(jié)果造成嚴(yán)重影響。因此，在數(shù)據(jù)分析的過程中，需要正確恰當(dāng)?shù)奶幚砣笔е?，使得研究更加有意義。

國外對缺失數(shù)據(jù)現(xiàn)象的問題很早就有了研究。自20世紀(jì)40年代起，有關(guān)學(xué)者便開始了對缺失數(shù)據(jù)問題的初步研究，并強(qiáng)調(diào)處理缺失數(shù)據(jù)問題的重要性。從缺失機(jī)制方面來看，Rubin（1976）［1］首次提出了關(guān)于數(shù)據(jù)缺失的缺失機(jī)制問題，并將其分為三種，隨機(jī)缺失（MAR）、完全隨機(jī)缺失（MCAR）和非隨機(jī)缺失（NMAR）。Sinna S K、Saha K 和 Wang S J（2014）［2］利用半?yún)?shù)方法處理非單調(diào)缺失數(shù)據(jù)，并提出了在不可忽視機(jī)制（NI）基礎(chǔ)下的NI-缺失機(jī)制。從處理缺失數(shù)據(jù)方法的方面來看，Rubin等人（1977）［3］首次提出了用來處理缺失數(shù)據(jù)的EM算法，這個算法的提出也是缺失數(shù)據(jù)處理方面一個極具意義的里程碑。

針對SAEM算法的研究，最初始于Marc Lavielle等人（1999）［4］，他們提出了一種新的方法，即隨機(jī)逼近EM（SAEM）算法，它通過一次隨機(jī)逼近過程的迭代代替了EM算法的期望步驟，并證明了在一些附加條件下，SAEM算法的吸引平穩(wěn)點對應(yīng)于函數(shù)的局部極大值。SAEM算法的提出大大推動了EM算法的發(fā)展。Wei Jiang等人（2020）［5］提出了EM算法的一個基于Metropolis-Hastings抽樣的隨機(jī)逼近版本（SAEM），對含有不完全數(shù)據(jù)的Logistic回歸進(jìn)行統(tǒng)計推斷。

從近十幾年來的發(fā)展來看，同其他數(shù)據(jù)缺失處理方法一樣，參數(shù)方法、非參數(shù)方法在缺失數(shù)據(jù)應(yīng)用方面發(fā)展更好。Yi-Hau Chen（2003）［6］等人提出了一種新的半?yún)?shù)估計器，使用加權(quán)的經(jīng)驗協(xié)變量分布（權(quán)重由回歸模型確定）來估計得分方程。Morikawa等人（2017）［7］提出了一種針對MNAR（非隨機(jī)缺失）數(shù)據(jù)的半?yún)?shù)極大似然方法，模型假設(shè)的響應(yīng)傾向部分使用了參數(shù)假設(shè)，而結(jié)果部分使用了非參數(shù)模型。龍兵，候蘭寶（2020）［8］針對定時截尾試驗的弊端提出了一個新的壽命試驗方案，基于試驗數(shù)據(jù)得到了似然函數(shù)，運用極大似然法得到了參數(shù)的迭代方程。

在眾多缺失機(jī)制中，NI-缺失機(jī)制與除缺失協(xié)變量自身之外的所有變量有關(guān)，這種缺失機(jī)制能夠充分的利用數(shù)據(jù)中含有的信息，由此來更好的處理缺失數(shù)據(jù)。SAEM算法具有效率高、效果好的優(yōu)點。本文通過引入NI-缺失機(jī)制來改進(jìn)SAEM算法，并通過模擬實驗，在同缺失率同缺失機(jī)制下，將其處理缺失數(shù)據(jù)的性能與半?yún)?shù)方法進(jìn)行對比，除標(biāo)準(zhǔn)誤差外，將回判準(zhǔn)確率作為評價兩種算法性能的標(biāo)準(zhǔn)。最后，將SAEM算法引入非酒精性脂肪肝數(shù)據(jù)并與半?yún)?shù)方法進(jìn)行對比。

1 SAEM算法

1.1 EM算法

SAEM的前身起源于1977年Rubin等人提出的EM算法，該算法中每一次迭代都分兩步，期望步（E步）和極大步（M步），從不完全數(shù)據(jù)中獲得極大似然估計量。給定初始值為θ0，xmis表示缺失的協(xié)變量，xobs表示觀察到的協(xié)變量，迭代k次后，通過以下兩步將θk-1更新為θk：

1.2 SAEM算法

1999年 Marc Lavielle，Bernard Delyon和 Eric Moulines首次提出隨機(jī)逼近EM（SAEM）算法，它通過一次隨機(jī)逼近過程的迭代代替了EM算法的期望步驟：

2 NI-缺失機(jī)制

定義示性函數(shù)為R，則有：

當(dāng)Rij=1的概率取決于完全觀察變量和部分觀察變量的值時，就會發(fā)生不可忽視（NI）機(jī)制，這意味著Xij的缺失機(jī)制可能取決于Xi=Xi1,…,Xip以及Yi和Zi的所有成分。對于NI缺失機(jī)制下的缺失數(shù)據(jù)，需要通過聯(lián)合數(shù)據(jù)生成過程和缺失過程的極大似然值來估計模型參數(shù)。NI缺失機(jī)制的缺點是數(shù)據(jù)會很少或根本沒有提供關(guān)于缺失模型中參數(shù)的信息。因此，模型參數(shù)是弱可辨識的。

為解決這一問題，Samiran Sinha等人（2014）［2］引入了另一種機(jī)制，用來處理可識別性問題。為處理缺失數(shù)據(jù)的單調(diào)模式，當(dāng)Rij=1的概率可能取決于除了第j個協(xié)變量Xij的所有變量時，定義缺失機(jī)制“NI-”，也就是說，P(Rij=1|Yi,Xi,Zi)=P(Rij=1|Yi,Xi(-j),Zi)。因此，NI-機(jī)制假設(shè)Rij取決于完全觀察到或部分缺失的其他變量，但不依賴于第j個協(xié)變量Xij的值，假設(shè)NI-缺失機(jī)制下，協(xié)變量具有以下的性質(zhì)：（1）Ri1,...,Rip獨立于Xi,Yi,Zi，（2）對于每個 j，P(Rij=1|Yi,Xi,Zi)>0對于每一個組合 Yi，Xi，Zi的概率為一。

根據(jù)NI-機(jī)制的特點可以得出以下結(jié)論：（1）NI-機(jī)制有助于減少不可忽視的缺失數(shù)據(jù)所遇到的可識別性問題。因此，該缺失機(jī)制下不需要對缺失協(xié)變量的分布進(jìn)行參數(shù)模型假設(shè)。（2）NI-機(jī)制允許估計模型參數(shù)，而無需為缺失協(xié)變量指定參數(shù)模型。當(dāng)假設(shè)為MAR缺失機(jī)制出現(xiàn)問題時，可以嘗試NI-缺失機(jī)制。（3）相比于其他缺失機(jī)制，NI-機(jī)制更能夠充分利用數(shù)據(jù)中已有的信息，以此來更加有效地對缺失數(shù)據(jù)進(jìn)行處理。

3 改進(jìn)的SAEM算法

2020年 Wei Jiang等人［5］的實驗中使用的是完全隨機(jī)缺失（MCAR）機(jī)制，MCAR缺失機(jī)制是指數(shù)據(jù)的缺失與所有變量都是無關(guān)的，只與常數(shù)項有關(guān)，這種缺失機(jī)制無法最大限度的利用觀測到的數(shù)據(jù)。Sinna等人提出的NI-缺失機(jī)制與除缺失協(xié)變量自身之外的所有變量有關(guān)，這種缺失機(jī)制能夠充分的利用數(shù)據(jù)中含有的信息。本文通過引入NI-機(jī)制來改進(jìn)SAEM算法，以此來更充分的利用數(shù)據(jù)中的信息，更好的處理缺失數(shù)據(jù)。改進(jìn)流程如下。

改進(jìn)的SAEM算法：

（1）NI-缺失機(jī)制：引入NI-缺失機(jī)制，在該缺失機(jī)制下生成缺失數(shù)據(jù)。

（2）模 擬 ：對 于 i=1,…,n，從 P(xi,mis|xi,obs；yi;θk-1)中得到。

（3）隨機(jī)逼近：根據(jù)下式更新函數(shù)Q：

其中（γk）是正數(shù)的非遞增序列。

（4）最大化：更新θ的估計：

4 模擬研究

4.1 基于NI-缺失機(jī)制下的SAEM算法

為了驗證NI-缺失機(jī)制的性能，本文通過NI-缺失機(jī)制生成缺失數(shù)據(jù)，并將生成的缺失數(shù)據(jù)作為樣本量，利用SAEM算法進(jìn)行回歸分析。

首先生成一組樣本量為n=200，協(xié)變量為p=5的模擬數(shù)據(jù)，然后根據(jù)NI-缺失機(jī)制來對數(shù)據(jù)進(jìn)行缺失，假設(shè)該模擬數(shù)據(jù)中所有協(xié)變量均等量缺失，且協(xié)變量之間不相關(guān)。設(shè)真正的參數(shù)值為：β=（-0.2，0.5，-0.3，1，0，-0.6）；五個協(xié)變量服從正態(tài)分布，其均值為：μ=（1，2，3，4，5），標(biāo)準(zhǔn)差為：σ=（1，2，3，4，5），設(shè)置 4組缺失率分別為10%、20%、30%和40%的數(shù)據(jù)。最后將不同缺失率下SAEM算法估計出來的參數(shù)值β的標(biāo)準(zhǔn)誤差STDβ作為判別標(biāo)準(zhǔn)之一，再利用估計得到的參數(shù)值β重新帶入Logistic回歸模型中，將得到的結(jié)果變量y值與y的真實值進(jìn)行對比，計算兩者的重復(fù)率，將其作為檢驗性能的第二個判別標(biāo)準(zhǔn)，并與2020年Wei Jiang等人使用的MCAR缺失機(jī)制進(jìn)行比較。兩種缺失機(jī)制下得到的標(biāo)準(zhǔn)誤差如表1所示，標(biāo)準(zhǔn)誤差表示為STDβ。

表1 兩種缺失機(jī)制下不同缺失率模擬結(jié)果

表1為不同缺失率下NI-缺失機(jī)制與MCAR缺失機(jī)制下參數(shù)估計結(jié)果的標(biāo)準(zhǔn)誤差，由表2可以看出，隨著缺失率的增長，NI-缺失機(jī)制下估計的標(biāo)準(zhǔn)誤差大多數(shù)情況下比MCAR缺失機(jī)制下估計的標(biāo)準(zhǔn)誤差更小。

表2為NI-缺失機(jī)制與MCAR缺失機(jī)制在不同缺失率下的回判結(jié)果，由表2可以得出結(jié)論：隨著缺失率的增長，NI-缺失機(jī)制下對y值進(jìn)行回判得到的準(zhǔn)確率比MCAR缺失機(jī)制下得到的更高。

表2 兩種缺失機(jī)制下不同缺失率的回判結(jié)果

除標(biāo)準(zhǔn)誤差與回判準(zhǔn)確率外，本文還使用了相對誤差對兩種缺失機(jī)制進(jìn)行了分析，分析結(jié)果與前兩種判別方法所得到的結(jié)果類似。不同缺失率下NI-缺失機(jī)制的相對誤差始終小于MCAR缺失機(jī)制下的相對誤差，MCAR缺失機(jī)制下相對誤差的極大值為92.36%，極小值為4.42%；NI-缺失機(jī)制下相對誤差的極大值為83.91%，極小值為2.44%。

4.2 NI-缺失機(jī)制下SAEM算法與半?yún)?shù)方法的對比

眾多缺失數(shù)據(jù)的處理方法中，半?yún)?shù)方法適用多個變量帶有缺失項的情況，與其他缺失數(shù)據(jù)處理方法相比較，半?yún)?shù)方法更能在不損失數(shù)據(jù)信息的前提下利用缺失部分的信息。近幾年來，由于半?yún)?shù)方法處理缺失數(shù)據(jù)的效果較好，半?yún)?shù)方法得到了較為廣泛的應(yīng)用。本節(jié)將半?yún)?shù)方法與引入NI-缺失機(jī)制的SAEM算法進(jìn)行對比分析。

首先生成一組樣本量為n=200，協(xié)變量為p=10的模擬數(shù)據(jù)，協(xié)變量分別服從分布：X1～Bernoulli(0.5)，當(dāng)r=2，3，4，5時，Xr～Normal(1，0.5)，當(dāng)r=6，7時 ，X～rNormal(2,2 )，X8～Normal(1,1)，X9=I(X2=1) Normal(- 0.25,1 ) +I(X2=0)Normal(0 .25,1 )，X10～Normal(-0.2X2+0.2X3，1)，

其中X1為完整數(shù)據(jù)，X2，…，X10為NI-缺失機(jī)制下等量缺失的不完全協(xié)變量，I是示性函數(shù)。生成四組缺失率分別為5%、10%、20%、30%和35%的數(shù)據(jù)，真正的參數(shù)值為：β=（-2.0，0.1，0.2，0.3，0.4，0.5，0.5，-0.5，-0.4，-0.3，-0.2）。將引入NI-缺失機(jī)制的SAEM算法和半?yún)?shù)方法分別對缺失數(shù)據(jù)進(jìn)行參數(shù)估計，并將估計出的參數(shù)值代入Logistic模型中，可以得到估計結(jié)果y值，對y值進(jìn)行回判所得到的準(zhǔn)確率如表3所示。

表3 不同缺失率下兩種方法的回判準(zhǔn)確率

由表3中的回判準(zhǔn)確率可以看出，當(dāng)缺失率為5%到20%之間時，SAEM算法比半?yún)?shù)方法估計結(jié)果的回判準(zhǔn)確率高；當(dāng)缺失率大于20%時，半?yún)?shù)的回判準(zhǔn)確率要更高。

不同缺失率下兩種方法模擬結(jié)果的標(biāo)準(zhǔn)誤差如表4、表5所示。從表4可以看出，當(dāng)缺失率為5%的情況下，SAEM算法得到的標(biāo)準(zhǔn)誤差比半?yún)?shù)得到的標(biāo)準(zhǔn)誤差要小；當(dāng)缺失率為10%時，兩者差距不大，約在0.01到0.05之間；從表5可以看出，當(dāng)缺失率為20%～30%的情況下，半?yún)?shù)方法的標(biāo)準(zhǔn)誤差小于SAEM算法，兩者差距小于約0.2；當(dāng)缺失率為35%時，半?yún)?shù)方法的標(biāo)準(zhǔn)誤差小于SAEM算法，兩者差距達(dá)到約0.3以上。由表4、表5可以得出結(jié)論：當(dāng)缺失率為5%時，SAEM算法對缺失數(shù)據(jù)參數(shù)估計的性能優(yōu)于半?yún)?shù)方法，當(dāng)缺失率大于10%時，半?yún)?shù)方法對缺失數(shù)據(jù)參數(shù)估計的性能優(yōu)于SAEM算法，但兩者差別不大。當(dāng)缺失率為35%時，兩者差別較大。

表4 缺失率小于20%的情況下兩種方法模擬結(jié)果的標(biāo)準(zhǔn)誤差

表5 缺失率≥20%的情況下兩種方法模擬結(jié)果的標(biāo)準(zhǔn)誤差

除標(biāo)準(zhǔn)誤差與回判準(zhǔn)確率外，本文還使用相對誤差對兩種算法進(jìn)行了比較分析，結(jié)果與表4，表5類似。當(dāng)缺失率小于10%時，SAEM算法的標(biāo)準(zhǔn)誤差范圍為0.77%～66.77%，半?yún)?shù)方法的標(biāo)準(zhǔn)誤差范圍是1.70%～85.45%，SAEM算法的標(biāo)準(zhǔn)誤差小于半?yún)?shù)方法；當(dāng)缺失率為10%～20%的情況下，SAEM算法的標(biāo)準(zhǔn)誤差范圍為6.24%～69.59%，半?yún)?shù)方法的標(biāo)準(zhǔn)誤差范圍是8.89%～58.37%，兩種算法的相對誤差差別不大，當(dāng)缺失率大于20%時，SAEM算法的標(biāo)準(zhǔn)誤差范圍為9.75%～81.03%，半?yún)?shù)方法的標(biāo)準(zhǔn)誤差范圍是9.73%～77.70%，半?yún)?shù)方法的相對誤差小于SAEM算法。

兩種算法在處理缺失數(shù)據(jù)時所需時間對比如表6所示。

表6 兩種算法運行時間對比

由表6可以看出，運行時間上SAEM算法遠(yuǎn)小于半?yún)?shù)方法，半?yún)?shù)方法的運行時間是SAEM算法的近669倍。

本節(jié)通過模擬實驗，從標(biāo)準(zhǔn)誤差、回判準(zhǔn)確率和運行時間方面將NI-缺失機(jī)制下的SAEM算法和半?yún)?shù)方法的性能進(jìn)行了對比分析。根據(jù)實驗數(shù)據(jù)可以得出結(jié)論：當(dāng)缺失率較小時，SAEM算法對缺失數(shù)據(jù)的處理性能優(yōu)于半?yún)?shù)方法，當(dāng)缺失率較大時，半?yún)?shù)方法的性能更好；在代碼運行過程中，SAEM算法的運行時間比半?yún)?shù)方法運行時間快600多倍，始終小于半?yún)?shù)方法。

5 實證分析

為了進(jìn)一步驗證兩種方法在處理不同缺失率下的缺失數(shù)據(jù)的性能，本節(jié)引入真實數(shù)據(jù)，將兩種方法分別進(jìn)行分析。

本節(jié)引入的真實數(shù)據(jù)來源于住院治療的患者，年齡均≥75歲，共116位患者作為本次研究的樣本量。樣本的結(jié)果變量即為非酒精性脂肪肝的患病狀態(tài)，設(shè)1表示患者患有非酒精性脂肪肝，0表示患者不患該疾病。樣本量中共有十個協(xié)變量，其中前七個協(xié)變量分別為：年齡、甘油三酯（TG）、總膽固醇（CHOL）、高密度脂蛋白膽固醇（HDL-C）、低密度脂蛋白膽固醇（LDL-C）、血清載脂蛋白A-1（ApoA1）、載脂蛋白B（ApoB），前七個變量為完全可觀測變量；后三個變量分別為：尿酸、胱抑素（Cys）、脂蛋白（a）（Lp（a）），后三個變量的缺失率分別為：20.69%，9.48%，13.79%。

將前七個可完全觀測到的變量表示為(Z1，Z2，Z3，Z4，Z5，Z6，Z7)，后 三 個 不 可 完 全觀測變量表示為：(X1，X2，X3)。根據(jù) Sinna.S 等人2014年［2］提出的文獻(xiàn)中鑒別NI-缺失機(jī)制的方法，設(shè)R1,R2,R3分別代表X1，X2，X3三個協(xié)變量的缺失指標(biāo)，當(dāng)Rx=1時，代表該協(xié)變量不缺失，當(dāng)Rx=0時，表示協(xié)變量缺失。在R1×R2=1的患者中，X2的缺失與X3密切相關(guān)（P值=0.004）。同樣，在R2×R3=1的患者中，X1的缺失與X3密切相關(guān)（P值=0.027），結(jié)果說明數(shù)據(jù)符合NI-缺失機(jī)制。根據(jù)兩種方法分別估計了參數(shù)，再將得到的參數(shù)值帶入Logistic模型，得到：

下面從回判的角度來比較兩種算法的性能。將兩種算法得到的回判準(zhǔn)確率的結(jié)果分別取十組數(shù)據(jù)，其中位數(shù)如表7所示。

表7 兩種算法對真實數(shù)據(jù)的回判準(zhǔn)確率

由表7可以看出，半?yún)?shù)方法的回判準(zhǔn)確率略大于SAEM算法?？赡艿脑蚴侨笔蔬^高（最高達(dá)到了20.69%），缺失率參差不齊，固在回判準(zhǔn)確率方面半?yún)?shù)優(yōu)于SAEM算法。但是從時間的角度來講，SAEM算法遠(yuǎn)優(yōu)于半?yún)?shù)方法。

6 結(jié)論

本文對SAEM算法進(jìn)行了改進(jìn)，引入了一種NI-缺失機(jī)制，以此來更充分的利用數(shù)據(jù)中的信息，更好的處理缺失數(shù)據(jù)。通過模擬實驗證明了該缺失機(jī)制下SAEM算法對缺失數(shù)據(jù)的處理性能比原始缺失機(jī)制的效果要更好，隨后將改進(jìn)的SAEM算法與半?yún)?shù)方法在模擬數(shù)據(jù)中進(jìn)行了比較，并將參數(shù)估計值、標(biāo)準(zhǔn)誤差、相對誤差和回判準(zhǔn)確率作為判別標(biāo)準(zhǔn)，最終得出結(jié)論：在低缺失率的數(shù)據(jù)中，SAEM算法對缺失數(shù)據(jù)的處理性能要優(yōu)于半?yún)?shù)方法，當(dāng)數(shù)據(jù)為高缺失率時，半?yún)?shù)方法要優(yōu)于SAEM算法。從效率方面來看，SAEM算法的運行時間要遠(yuǎn)小于半?yún)?shù)方法。在大數(shù)據(jù)時代，在線對具有缺失值的數(shù)據(jù)進(jìn)行研究，效率顯得尤為主要。因此，本文所提出的NI-缺失機(jī)制下的SAEM算法更具有優(yōu)勢。

本文所做的對比研究都是基于正態(tài)分布下的協(xié)變量缺失的情況，在非正態(tài)分布情況下的結(jié)果以及令SAEM算法與似然函數(shù)相結(jié)合等方面值得未來進(jìn)一步研究。