龐新生

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

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

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

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

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

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

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

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

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

定義函數(shù):

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

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

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

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

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

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

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

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

2.1 MCMC方法

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

2.1.1 插補(bǔ)步驟

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

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

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

2.1.2 后驗(yàn)步驟

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

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

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

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

2.2 DA算法

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

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

圖1 DA算法迭代模擬過(guò)程

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

圖2 快速收斂(時(shí)序圖)

圖3 緩慢收斂(時(shí)序圖)

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

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

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

[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)計(jì)決策論及貝葉斯分析[M].北京：中國(guó)統(tǒng)計(jì)出版社，1997.

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