基于SAEM算法對缺失協變量的Logistic模型參數估計

劉玥，施三支

（長春理工大學 理學院，長春 130022）

數據的缺失幾乎存在于所有的實證研究領域。可能會出現數據缺失的原因有很多，包括調查無響應、無法測量和單純的數據丟失等。但同時，現在各類研究不斷増多，這就需要不斷采集挖掘數據；數據量越大，數據缺失的可能性也就越大，缺失數據就會對統計推斷結果造成嚴重影響。因此，在數據分析的過程中，需要正確恰當的處理缺失值，使得研究更加有意義。

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

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

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

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

1 SAEM算法

1.1 EM算法

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

1.2 SAEM算法

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

2 NI-缺失機制

定義示性函數為R，則有：

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

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

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

3 改進的SAEM算法

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

改進的SAEM算法：

（1）NI-缺失機制：引入NI-缺失機制，在該缺失機制下生成缺失數據。

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

（3）隨機逼近：根據下式更新函數Q：

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

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

4 模擬研究

4.1 基于NI-缺失機制下的SAEM算法

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

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

表1 兩種缺失機制下不同缺失率模擬結果

表1為不同缺失率下NI-缺失機制與MCAR缺失機制下參數估計結果的標準誤差，由表2可以看出，隨著缺失率的增長，NI-缺失機制下估計的標準誤差大多數情況下比MCAR缺失機制下估計的標準誤差更小。

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

表2 兩種缺失機制下不同缺失率的回判結果

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

4.2 NI-缺失機制下SAEM算法與半參數方法的對比

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

首先生成一組樣本量為n=200，協變量為p=10的模擬數據，協變量分別服從分布：X1～Bernoulli(0.5)，當r=2，3，4，5時，Xr～Normal(1，0.5)，當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為完整數據，X2，…，X10為NI-缺失機制下等量缺失的不完全協變量，I是示性函數。生成四組缺失率分別為5%、10%、20%、30%和35%的數據，真正的參數值為：β=（-2.0，0.1，0.2，0.3，0.4，0.5，0.5，-0.5，-0.4，-0.3，-0.2）。將引入NI-缺失機制的SAEM算法和半參數方法分別對缺失數據進行參數估計，并將估計出的參數值代入Logistic模型中，可以得到估計結果y值，對y值進行回判所得到的準確率如表3所示。

表3 不同缺失率下兩種方法的回判準確率

由表3中的回判準確率可以看出，當缺失率為5%到20%之間時，SAEM算法比半參數方法估計結果的回判準確率高；當缺失率大于20%時，半參數的回判準確率要更高。

不同缺失率下兩種方法模擬結果的標準誤差如表4、表5所示。從表4可以看出，當缺失率為5%的情況下，SAEM算法得到的標準誤差比半參數得到的標準誤差要小；當缺失率為10%時，兩者差距不大，約在0.01到0.05之間；從表5可以看出，當缺失率為20%～30%的情況下，半參數方法的標準誤差小于SAEM算法，兩者差距小于約0.2；當缺失率為35%時，半參數方法的標準誤差小于SAEM算法，兩者差距達到約0.3以上。由表4、表5可以得出結論：當缺失率為5%時，SAEM算法對缺失數據參數估計的性能優于半參數方法，當缺失率大于10%時，半參數方法對缺失數據參數估計的性能優于SAEM算法，但兩者差別不大。當缺失率為35%時，兩者差別較大。

表4 缺失率小于20%的情況下兩種方法模擬結果的標準誤差

表5 缺失率≥20%的情況下兩種方法模擬結果的標準誤差

除標準誤差與回判準確率外，本文還使用相對誤差對兩種算法進行了比較分析，結果與表4，表5類似。當缺失率小于10%時，SAEM算法的標準誤差范圍為0.77%～66.77%，半參數方法的標準誤差范圍是1.70%～85.45%，SAEM算法的標準誤差小于半參數方法；當缺失率為10%～20%的情況下，SAEM算法的標準誤差范圍為6.24%～69.59%，半參數方法的標準誤差范圍是8.89%～58.37%，兩種算法的相對誤差差別不大，當缺失率大于20%時，SAEM算法的標準誤差范圍為9.75%～81.03%，半參數方法的標準誤差范圍是9.73%～77.70%，半參數方法的相對誤差小于SAEM算法。

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

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

由表6可以看出，運行時間上SAEM算法遠小于半參數方法，半參數方法的運行時間是SAEM算法的近669倍。

本節通過模擬實驗，從標準誤差、回判準確率和運行時間方面將NI-缺失機制下的SAEM算法和半參數方法的性能進行了對比分析。根據實驗數據可以得出結論：當缺失率較小時，SAEM算法對缺失數據的處理性能優于半參數方法，當缺失率較大時，半參數方法的性能更好；在代碼運行過程中，SAEM算法的運行時間比半參數方法運行時間快600多倍，始終小于半參數方法。

5 實證分析

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

本節引入的真實數據來源于住院治療的患者，年齡均≥75歲，共116位患者作為本次研究的樣本量。樣本的結果變量即為非酒精性脂肪肝的患病狀態，設1表示患者患有非酒精性脂肪肝，0表示患者不患該疾病。樣本量中共有十個協變量，其中前七個協變量分別為：年齡、甘油三酯（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)。根據 Sinna.S 等人2014年［2］提出的文獻中鑒別NI-缺失機制的方法，設R1,R2,R3分別代表X1，X2，X3三個協變量的缺失指標，當Rx=1時，代表該協變量不缺失，當Rx=0時，表示協變量缺失。在R1×R2=1的患者中，X2的缺失與X3密切相關（P值=0.004）。同樣，在R2×R3=1的患者中，X1的缺失與X3密切相關（P值=0.027），結果說明數據符合NI-缺失機制。根據兩種方法分別估計了參數，再將得到的參數值帶入Logistic模型，得到：

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

表7 兩種算法對真實數據的回判準確率

由表7可以看出，半參數方法的回判準確率略大于SAEM算法。可能的原因是缺失率過高（最高達到了20.69%），缺失率參差不齊，固在回判準確率方面半參數優于SAEM算法。但是從時間的角度來講，SAEM算法遠優于半參數方法。

6 結論

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

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