缺失數(shù)據(jù)多重插補(bǔ)處理方法的算法實現(xiàn)

龐新生

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

1 EM方法在缺失數(shù)據(jù)多重插補(bǔ)處理中的實現(xiàn)

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

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

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

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

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

其中θ()i為已知的當(dāng)前參數(shù)的估計值。

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

定義函數(shù):

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

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

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

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

其中Θ代表參數(shù)空間。

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

2 MCMC方法在缺失數(shù)據(jù)多重插補(bǔ)處理中的實現(xiàn)

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

2.1 MCMC方法

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

2.1.1 插補(bǔ)步驟

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

其中∑11是Xobs的協(xié)方差矩陣，∑22是Xmis的協(xié)方差矩陣，∑12是Xobs與Xmis之間的協(xié)方差矩陣。在多元正態(tài)分布的假設(shè)下，當(dāng)給定Xobs=x1時，Xmis的均數(shù)為：

其對應(yīng)的條件協(xié)方差矩陣為：

2.1.2 后驗步驟

在每一次循環(huán)運(yùn)算中，用上一次插補(bǔ)步中得到的μ和協(xié)方差矩陣對參數(shù)φ進(jìn)行模擬。循環(huán)進(jìn)行這兩步過程，產(chǎn)生一個足夠長的馬爾科夫鏈：

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

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

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

2.2 DA算法

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

使用DA去實現(xiàn)多重插補(bǔ)，為了產(chǎn)生恰當(dāng)?shù)亩嘀夭逖a(bǔ)，我們必須從數(shù)據(jù)增廣中迭代Ymis形成或者形成m條長度為t獨(dú)立鏈。為了估計的需要，必須有參數(shù)的初始值，通過EM進(jìn)行ML估計的結(jié)果是一個很好的選擇。同時應(yīng)該注意的是，必須需要選擇一個比較大的t以確保連續(xù)插補(bǔ)統(tǒng)計上的獨(dú)立。

圖1 DA算法迭代模擬過程

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

圖2 快速收斂(時序圖)

圖3 緩慢收斂(時序圖)

圖4 快速收斂(自相關(guān)圖)

圖5 緩慢收斂(自相關(guān)圖)

從MCMC方法(或DA算法)的思想可以看出，基于模擬思想的多重插補(bǔ)也可以用于處理單位無回答，此時，只需要模擬含缺失數(shù)據(jù)變量或參數(shù)的聯(lián)合分布，進(jìn)行隨機(jī)取值，從而創(chuàng)建插補(bǔ)數(shù)據(jù)集。各個插補(bǔ)數(shù)據(jù)集分析結(jié)果的合并也遵循多重插補(bǔ)推論和Rubin的合并規(guī)則。作為計算方法，MCMC方法(或DA算法)也存在一些不足之處：一是需要多元正態(tài)假設(shè)；二是過程復(fù)雜運(yùn)算繁瑣；三是對于是否收斂不好確定。慶幸的是SAS、S-PLUS、MICE中提供了MCMC運(yùn)算，使得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著，賈乃光譯.統(tǒng)計決策論及貝葉斯分析[M].北京：中國統(tǒng)計出版社，1997.

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