缺失數據多重插補處理方法的算法實現

龐新生

在國外相當多的抽樣調查中，對缺失值進行插補處理是非常普遍的，替換缺失數據技術的意義在于比列表刪除浪費更少的信息，且當缺失數據為非隨機缺失時，替換缺失數據技術比列表刪除更穩健，特別是當數據收集者與數據分析者是不同的個體時，插補法更具優勢。插補法主要經歷了單一插補和多重插補兩階段，多重插補法的出現，彌補了單一插補法的缺陷，第一，多重插補過程產生多個中間插補值，可以利用插補值之間的變異反映無回答的不確定性，包括無回答原因已知情況下抽樣的變異性和無回答原因不確定造成的變異性。第二，多重插補通過模擬缺失數據的分布，較好地保持變量之間的關系。第三，多重插補能給出衡量估計結果不確定性的大量信息，單一插補給出的估計結果則較為簡單。與單一插補相比，多重插補唯一的缺點是需要做大量的工作來創建插補集并進行結果分析，無論是何種情況下的多重插補，其處理過程都是比較復雜的，新的統計計算方法的出現大大簡化了計算并完成一系列簡單的極大化或模擬。在缺失數據處理中，主要涉及的是數據添加算法，其中討論最多的是EM算法和馬爾科夫鏈蒙特卡洛方法(MCMC)。

1 EM方法在缺失數據多重插補處理中的實現

EM算法是Dempster，Laired和Rubin于1977年提出的求參數極大似然估計或最大后驗估計的一種方法，通過假設潛在變量的存在，EM算法極大地簡化了似然函數，從而解決了方程求解問題。

假設X是服從某一分布的觀測數據集，Y為缺失數據，則有完全數據集Z=(X,Y)，則Z的密度函數為:

從式(1)可以看出，密度函數 p(z|θ)是由邊緣密度函數p(x|θ)，缺失數據 y的假設，參數θ初始估計值及缺失數據與觀測變量之間的關系決定的。由式(11)給出的密度函數可以定義完全數據似然函數:

由于缺失數據未知，因此似然函數L(θ|Z)是隨機的，且由缺失數據Y所決定的。這里，我們假定存在缺失數據的變量 是隨機缺失的（MAR），在此假定之下，可以保證似然估計的精度。

由于E步是在給定觀測X和當前參數估計值，計算完全數據對數似然函數log p(X,Y|θ)關于缺失數據Y的期望，為此，定義對數似然函數的期望:

其中θ()i為已知的當前參數的估計值。

在式(3)中，X與θ()i為常數，θ為待優化的參數，Y為一隨機變量，并假設它服從某一分布 fY()y:

定義函數:

h()θ,Y=Δln p(X,Y|θ)(5)因此，式(5)可寫為:

其中f(y | X，θ(i))是缺失數據Y的邊緣密度函數，并且依賴于觀測數據和當前參數θ()i，D為 y的取值空間。由于有:

且因子 f(X |θ(i))與θ無關，所以在實際問題處理中，用f(y , X|θ(i))代替f(y | X，θ(i))不影響式(5)中似然函數的最優化。

EM算法的第二步M-step：最大化期望值 g(θ|θ(i),X)，即找到一個θ(i+1)，滿足：

其中Θ代表參數空間。

EM算法是利用缺失數據和模型參數之間的迭代關系：如果缺失數據已知，模型參數未知，那么根據缺失數據就可以對模型參數進行估計。同樣，如果模型參數已知，根據模型也可以得到缺失數據的估計。先在假定模型參數的基礎上得到缺失數據的估計，然后再利用缺失值的估計值修正模型參數，這樣不斷地進行迭代，直至模型參數值收斂。EM算法的主要目的在于提供一個簡單的迭代算法來計算極大似然估計，每一步迭代都能保證似然函數值增加，并且收斂到一個局部極大值，該算法的最大優點是簡單和穩定，不足之處在于：第一，不同的模型需要不同的程序，不存在統一的處理程序；第二，當缺失數據比較多時，運算速度往往比較慢；第三，標準差只能在運算收斂后通過其他均值計算，無法直接獲得。

2 MCMC方法在缺失數據多重插補處理中的實現

MCMC方法適合于處理多維非單調確定缺失目標變量多重插補，實踐中，一般都是通過DA法實現對復雜分布的模擬，從而使得多重插補得以實施。MCMC方法是一組方法的集合，最早用于探測分子布朗運動的分布。MCMC方法是通過運行很長一段時間后形成Markov鏈樣本，以便用樣本均值近似地求數學期望。構造這種Markov鏈的方法較多，其中包括Gibbs抽樣在內，大都是Metroplis-Hasting算法的特例，MCMC方法實質上就是利用Markov鏈進行Monte Carlo積分，在利用通用軟件來分析許多復雜的問題時，MCMC方法可提供統一的結構框架，在多重插補中旨在通過馬爾科夫鏈使缺失數據和參數的分布收斂，從而模擬其分布。

2.1 MCMC方法

MCMC是貝葉斯推斷中一種探索后驗分布的方法，下面通過正態模型說明MCMC基本思想和實施步驟，令Y=(y1,y2,…,yn)T為完全數據集，假定 y1，y2，…，yn|θ～iid Np(μ ，∑ )，其中 θ=(μ ，∑ )未知，運用該方法對該缺失數據集插補可以分為兩步：

2.1.1 插補步驟

根據給定的均數向量μ和協方差矩陣∑，從條件分布 p(Xmis|Xobs,φ)中為缺失值抽取插補值。假設是兩部分變量的均數向量，μ1是Xobs的均值向量，μ2是Xmis的均值向量，同時設定：

其中∑11是Xobs的協方差矩陣，∑22是Xmis的協方差矩陣，∑12是Xobs與Xmis之間的協方差矩陣。在多元正態分布的假設下，當給定Xobs=x1時，Xmis的均數為：

其對應的條件協方差矩陣為：

2.1.2 后驗步驟

在每一次循環運算中，用上一次插補步中得到的μ和協方差矩陣對參數φ進行模擬。循環進行這兩步過程，產生一個足夠長的馬爾科夫鏈：

當該鏈會聚在一個穩定的分布 p(Xmis,φ|Xobs)時，就可以近似獨立地從該分布中為缺失值抽取插補值。

為了建立插補程序，我們必須做某些假定：首先要求對缺失機制必須做出假定，如隨機缺失（MAR），如同可忽略的假定，令Yobs為觀測值，Ymis為缺失值，R為回答指示變量，R的分布依賴于Yobs而不依賴于Ymis，即有P(R |Yobs，Ymis)=P(R |Yobs)；其次要求對參數的先驗分布必須做出假定，多重插補必須反映插補模型參數的不確定性：

其中有：P()θ|Yobs∝L()θ|Yobsπ()θ，對于先驗分布π()θ要求，推斷對于π的選擇不敏感。

2.2 DA算法

MCMC方法構造馬氏鏈去模擬后驗分布f(Ymis|Yobs)，可以通過DA算法實現，該方法是MCMC算法之一，特別適合于缺失數據的處理。DA算法的特點在于可以處理任意缺失模式，具體插補過程如圖1所示。DA算法經過t次迭代后收斂于一個分布而不是一個值，收斂速度與數據缺失程度相關，如果數據缺失嚴重，收斂速度很慢，迭代的次數要多些，反之，收斂速度很快。總的來說，DA算法是重復兩個步驟，即：I步，從Pr(Ymis|Yobs,θ(t))中抽取中抽取θ(t+1)。重復該過程多次，這樣就建立了一條markovchain而該鏈收斂于P(Ymis，θ|Yobs)，這個分布就是對實際分布的預測。DA法估計的目的是從收斂的分布中隨機抽取Ymis值，替代缺失數據。當有關于均值向量和協方差矩陣的先驗信息時，直接利用先驗信息，就可以進行迭代。當先驗信息缺失時，利用大樣本理論，可以認為協方差矩陣∑服從∑(t+1)|Y～W-1(n-1,(n-1)S)的分布。均值向量矩陣U服從，其中W表示Wishart分布。

使用DA去實現多重插補，為了產生恰當的多重插補，我們必須從數據增廣中迭代Ymis形成或者形成m條長度為t獨立鏈。為了估計的需要，必須有參數的初始值，通過EM進行ML估計的結果是一個很好的選擇。同時應該注意的是，必須需要選擇一個比較大的t以確保連續插補統計上的獨立。

圖1 DA算法迭代模擬過程

運用DA算法時，為使各插補值盡量保持獨立，一般需迭代m×t次，得到Y(t)mis，Y(2t)mis，…，Y(mt)mis，這就是最終的m個插補值，其中t可以通過參數的時間序列圖和自相關函數圖(ACF)來確定，下面通過例子對這兩種方法分別說明。方法一，畫出參數θ與迭代次數的分布圖，即θ的時序圖，看其在何時趨于收斂，如果參數θ的時序圖沒有長期趨勢，我們稱之為快速收斂，如圖2所示，如果存在長期趨勢和變化，我們稱為緩慢收斂，如圖3所示；方法二，畫出參數θ的自相關圖(ACF)，自相關函數圖估計了每次迭代參數與k次迭代參數之間的相關系數，這些圖可以幫助數據分析人員去判斷經過多少次迭代后參數值與初始值之間就相互獨立了。每一個自相關函數圖顯示了一系列上下限值，在圖4、圖5上用兩條橫線表示，如果超出橫線，說明自相關系數是顯著的(α=0.05)。根據自相關系數收斂時的迭代次數，如果ACF很快衰減至0，我們稱之為快速收斂，如圖4所示，經過7次迭代后，ACF很快衰減至0；如果衰減很慢，我們稱之為緩慢收斂，如圖5所示，經過100次迭代后，ACF還沒有衰減至0。為了得到ACF的平穩估計，特別是當緩慢收斂時，需要相當多次迭代，而且從保守的角度來看，循環次數應該足夠大。一般情況下，希望自相關的初始值是正值，經過迭代很快或逐漸降為0，即使后面仍在迭代，其值仍然為0。為了提高收斂速度，在實施DA法之前，最好是先進行EM法的運算，DA算法通常以EM算法得到的結果作為初始值進行迭代。關于DA算法與EM算法之間的關系，有關研究給出了相應準則：如果EM算法經過t次迭代收斂，那么DA算法經過t次迭代幾乎也確定收斂。需要注意的是EM算法收斂于一個參數估計值，而DA算法收斂于參數值的分布。

圖2 快速收斂(時序圖)

圖3 緩慢收斂(時序圖)

圖4 快速收斂(自相關圖)

圖5 緩慢收斂(自相關圖)

從MCMC方法(或DA算法)的思想可以看出，基于模擬思想的多重插補也可以用于處理單位無回答，此時，只需要模擬含缺失數據變量或參數的聯合分布，進行隨機取值，從而創建插補數據集。各個插補數據集分析結果的合并也遵循多重插補推論和Rubin的合并規則。作為計算方法，MCMC方法(或DA算法)也存在一些不足之處：一是需要多元正態假設；二是過程復雜運算繁瑣；三是對于是否收斂不好確定。慶幸的是SAS、S-PLUS、MICE中提供了MCMC運算，使得MCMC越來越成為一種主流方法。

[1][美]Roderick J.A.Little,Donald B.Rubin Statistical Analysis with Missing Data[M].New York:John Wiley&Sons INC,2002.

[2][美]James O．Berger著，賈乃光譯.統計決策論及貝葉斯分析[M].北京：中國統計出版社，1997.

[3][美]Donald.B.Rubin Multiple Imputation For Nonresponse in Surveys[M].New Yrok:Jghn Wiley&Sons INC,1987.