于力超

（中央民族大學(xué) 理學(xué)院，北京 100081）

0 引言

調(diào)查過程中常出現(xiàn)數(shù)據(jù)缺失，在這種情形下不能直接將含缺失值的個體信息全部刪除而只利用數(shù)據(jù)完全的樣本數(shù)據(jù)，因為數(shù)據(jù)缺失與否往往與該數(shù)據(jù)的具體值系統(tǒng)地相關(guān)，比如關(guān)心的因變量是受訪者的年收入，高收入群體可能不愿意將自己的收入告知他人，這種情況造成的數(shù)據(jù)缺失就不能采用直接刪除法，否則會造成估計結(jié)果偏低。為此收集盡量多的協(xié)變量數(shù)據(jù)就很有意義?？梢圆捎貌逖a法（包括回歸插補法、近鄰插補法）或似然法（通過建立因變量、協(xié)變量和缺失指示變量的分布模型）[1]。在這些方法中，以往研究都假設(shè)協(xié)變量數(shù)據(jù)沒有缺失，回歸插補法通過有回答的樣本數(shù)據(jù)建立模型，代入因變量數(shù)據(jù)缺失樣本的協(xié)變量數(shù)據(jù)以估計因變量值。近鄰插補法通過定義距離函數(shù)，考慮用與因變量缺失個體的協(xié)變量最相似有數(shù)據(jù)的個體因變量值插補。似然法則引入數(shù)據(jù)缺失機制，通過建模得到因變量與缺失指示變量的聯(lián)合分布，進(jìn)而得到參數(shù)的極大似然估計。

實際調(diào)查中，協(xié)變量缺失的情況也是非常常見的，例如出于隱私考慮不愿公布年齡或婚姻狀況，由于知識水平所限無法回答某些超出其知識范圍的問題，也可能某些協(xié)變量收集難度大成本高，調(diào)查人員選擇只收集一部分受訪者的信息。協(xié)變量缺失時因變量可能沒有缺失，也可能協(xié)變量和因變量同時缺失，這就大大增加了協(xié)變量關(guān)于因變量回歸系數(shù)估計的難度。但目前國內(nèi)外關(guān)于協(xié)變量缺失情形下尤其是缺失機制是非隨機缺失時參數(shù)的估計方法研究還不多，Little（1992）[2]研究了回歸分析中缺失協(xié)變量的處理方法，但其研究局限于因變量與協(xié)變量之間滿足線性回歸關(guān)系，且假設(shè)因變量Y與協(xié)變量向量X的聯(lián)合分布為多元正態(tài)分布，Horton和Laird（1999）[3]研究了缺失機制為MAR的屬性協(xié)變量缺失問題，通過建立廣義線性模型，用極大似然法進(jìn)行參數(shù)估計，Michiels等（1999）[4]用選擇模型和模式混合模型建立似然函數(shù)，用極大似然法估計屬性協(xié)變量數(shù)據(jù)缺失時的總體參數(shù)。在軟件功能日益強大的今天，數(shù)據(jù)擴(kuò)充算法、EM算法、Gibbs抽樣法等迭代算法應(yīng)用日益廣泛。SAS Proc MI可以采用多重插補法進(jìn)行協(xié)變量數(shù)據(jù)缺失情形下的參數(shù)估計。R有多個軟件包可以進(jìn)行缺失數(shù)據(jù)統(tǒng)計分析，ACD包可以在因變量數(shù)據(jù)缺失時進(jìn)行屬性數(shù)據(jù)分析，mvnmle包在因變量和協(xié)變量聯(lián)合分布為多元正態(tài)分布時，進(jìn)行協(xié)變量數(shù)據(jù)缺失情形下的參數(shù)極大似然估計，MICE包是R中目前最常用的用于缺失數(shù)據(jù)分析的軟件包，可以進(jìn)行多變量缺失數(shù)據(jù)的多重插補，在多個協(xié)變量都可能存在缺失值時，使用MICE包中的mice函數(shù)，通過變量之間的關(guān)系預(yù)測缺失數(shù)據(jù)，利用蒙特卡洛方法生成多個完整數(shù)據(jù)集存在imp中，再對imp進(jìn)行線性回歸，最后用pool函數(shù)對回歸結(jié)果進(jìn)行匯總。

本文嘗試引入以上統(tǒng)計計算方法，采用多重插補法、Bayes法和極大似然法，在缺失機制為MAR或NMAR情形下，估計調(diào)查數(shù)據(jù)中協(xié)變量數(shù)據(jù)缺失情形下，協(xié)變量關(guān)于因變量影響大小的回歸系數(shù)。

1 符號、概念及協(xié)變量缺失處理方法

1.1 符號表示

設(shè)共有N個受訪者，因變量為Y，得到觀測值yi，i=1，...，N,其中可能有缺失值也可能沒有，同時對第i個受訪者，調(diào)查與yi相關(guān)的p維協(xié)變量Xi,Xi=(xi1，...，xip),Nxp 維協(xié)變量矩陣X=(X1，...，XN)'，假設(shè) (Xi，yi)，i=1,…,N獨立同分布。協(xié)變量含有缺失值的情形下，記第i個受訪者協(xié)變量的所有觀測值為Xobs，i，缺失值為Xmis，i，Xi=(Xobs，i，Xmis，i)， ，記 Nxp 維協(xié)變量缺失指示變量矩陣r=(r1，...，rN)',若第i個受訪者的第k個協(xié)變量數(shù)據(jù)未缺失rik=1,否則rik=0。

1.2 數(shù)據(jù)缺失機制的概念

本文考慮協(xié)變量數(shù)據(jù)缺失的情形，為表述簡便，假設(shè)收集到的樣本數(shù)據(jù)中，因變量Y無缺失，協(xié)變量X有缺失（若Y也有缺失，則只需在迭代過程中增加一步：從Ymis，i關(guān)于(yobs，i，Xi)后驗分布中抽取下一輪迭代的初值。除此沒有本質(zhì)上的差異）。此時定義的缺失機制與描述數(shù)據(jù)缺失與否與因變量數(shù)據(jù)之間的關(guān)系的常規(guī)定義不同，對于縱向缺失數(shù)據(jù)，缺失數(shù)據(jù)機制可分為完全隨機缺失（MCAR）、隨機缺失（MAR）和非隨機缺失（NMAR），MCAR指變量數(shù)據(jù)缺失與否與變量的值無關(guān)，這時候缺失數(shù)據(jù)集是總體數(shù)據(jù)集的一個簡單隨機樣本，僅用觀測數(shù)據(jù)就可得總體參數(shù)的無偏估計，協(xié)變量缺失時MAR指協(xié)變量是否缺失僅與已經(jīng)觀測到的協(xié)變量值有關(guān)，即L，可見在給定Xobs，i時，Xmis，i與ri獨立，即，可以忽略缺失值，僅用完全觀測數(shù)據(jù)對總體參數(shù)進(jìn)行估計，這就是MAR被稱為可忽略缺失機制的原因。

2 插補法處理缺失協(xié)變量

2.1 完全數(shù)據(jù)法

該方法將所有含缺失值的受訪者刪除，僅保留完全數(shù)據(jù)樣本進(jìn)行統(tǒng)計分析，這種方法僅當(dāng)缺失機制為MCAR時可以得到參數(shù)無偏估計，缺失機制為MAR或NMAR時得到的估計量有偏，若協(xié)變量數(shù)目較多，即使某受訪者只有一個協(xié)變量數(shù)據(jù)缺失，仍需要將其所有數(shù)據(jù)刪去，這樣可能導(dǎo)致可用的數(shù)據(jù)量不足，造成大量的信息浪費。雖然這種方法經(jīng)常采用，但并不值得推薦。

2.2 可用信息法

該方法利用所有可利用的信息進(jìn)行參數(shù)估計，若N個受訪者中，第k個協(xié)變量有觀測值的受訪者個數(shù)為Nk，則第k個協(xié)變量的期望μk就用這Nk個數(shù)據(jù)的均值估計；若N個受訪者中，第j個和第k個協(xié)變量都有觀測值的受訪者個數(shù)為Njk，則第j個協(xié)變量和第k個協(xié)變量之間的協(xié)方差σjk用這Njk個受訪者的相應(yīng)觀測數(shù)據(jù)計算。但此時關(guān)于第

由上可見，如何獲得協(xié)變量缺失值的插補值是關(guān)鍵，這就要求獲得后驗分布：

具體地，若協(xié)變量X均為離散型，則：

分母中是對第i個受訪者所有缺失協(xié)變量的所有可能取值求和。若協(xié)變量均為連續(xù)型，則分母中求和號改為積分號。一般情況下，協(xié)變量既有離散型又有連續(xù)型，記第i個受訪者離散型協(xié)變量缺失值的所有可能取值集合為，連續(xù)型含缺失值協(xié)變量集合記為，缺失值的可能取值集合為Ω，則：

由式（3）可見，為獲得插補值Xmis，i，需要獲得Xmis后驗分布P(Xmis，i|Xobs，i，yi，γ)的參數(shù)γ，然后從該后驗分布中抽取Xmis，i，獲得參數(shù)γ的方法有兩種:在軟件SAS Proc MI中，假設(shè)（Yi，Xobs，i，Xmis，i）服從多元正態(tài)分布，并假設(shè)協(xié)變量的缺失機制為MAR，由于MAR為可忽略缺失機制，可以僅通過觀測的數(shù)據(jù)對總體參數(shù)進(jìn)行估計，通過軟件可以得到P(Xmis，i|Xobs，i，yi，γ)中參數(shù)γ的極大似然估計?，代入Xmis，i的后驗分布，然后從P(Xmis，i|Xobs，i，yi，?)中多次抽取Xmis，i，從而得到多個插補數(shù)據(jù)集進(jìn)行參數(shù)估計。但這種方法需要做總體服從多元正態(tài)分布的假定。i個受訪者第j個協(xié)變量和第k個協(xié)變量之間的p階相關(guān)矩陣不一定正定，若非正定，協(xié)變量X與因變量Y之間的回歸系數(shù)β將無法確定。

2.3 多重插補法

采用多重插補法進(jìn)行參數(shù)估計的基本步驟如下[5]：第一步：抽取符合缺失變量后驗分布的插補值，重復(fù)進(jìn)行M次得到M個“完整”數(shù)據(jù)集；

第二步：對第m個插補后的“完整”數(shù)據(jù)集進(jìn)行參數(shù)估計，記待估參數(shù)為γ，得到估計值γ?(m)，其中m=1,…,M。第三步：參數(shù)估計為：

下面采用Gibbs抽樣獲得插補值的方法更加合理。

根據(jù)式（3），后驗分布P(Xmis，i|Xobs，i，yi，γ)中的參數(shù)γ包括P(yi|Xi，β)中的參數(shù)β和Xi邊緣分布中的參數(shù)α，Gibbs抽樣分兩步進(jìn)行:

第一步：從觀測值（Yi，Xobs，i）已知時（β，α）的后驗分布中抽?。é?，α）參數(shù)值，由于已知觀測值時（β，α）的后驗分布形式未知，實際操作中簡便起見常取無信息分布（例如定義域上的均勻分布）作為(β,α)的后驗分布。為了充分利用觀測值信息，建議通過數(shù)據(jù)擴(kuò)充算法(Data augmentation)[6]來實現(xiàn)對 (β,α)的 Gibbs抽樣。數(shù)據(jù)擴(kuò)充算法基本思想體現(xiàn)在下面兩個式子中：

其中S(Xmis，i|Xobs，i，Yi)和S(γ|Xobs，i，Yi)分別為在給定觀測值條件下，協(xié)變量缺失值和參數(shù)分布的支撐集。綜合以上兩式，可以得到參數(shù)后驗分布的迭代算法，設(shè)第k步迭代得到參數(shù)γ的后驗分布為fk(γk|Xobs，i，Yi)，則第 k+1步參數(shù)γ的后驗分布為：

Tanner和Wong(1987)證明了在一定正則條件下，fk+1依分布收斂于參數(shù)的后驗分布

根據(jù)以上兩式，本文建立了迭代算法：

首先，設(shè)定參數(shù)的初始值γ0開始迭代，在第k步參數(shù)后驗分布中抽取 m個參數(shù)值，對每

其次，將第k步的m個插補值代入數(shù)據(jù)集得到m個完整數(shù)據(jù)集（當(dāng)缺失機制為不可忽略時，無法僅通過觀測值（Yi，Xobs，i）得到總體參數(shù)的后驗分布，插補后通過完整數(shù)據(jù)集則可以得到參數(shù)的后驗分布，這是數(shù)據(jù)擴(kuò)充算法的意義所在），通過這m個完整數(shù)據(jù)集下參數(shù)的后驗分布可以得到第k+1步參數(shù)的后驗分布：

如此迭代進(jìn)行，可以得到參數(shù)γ的后驗分布中 取γ（r）代 入中 抽 取，如此重復(fù)進(jìn)行m次，得到m個完整數(shù)據(jù)集，進(jìn)而用對完整數(shù)據(jù)集用式（1）和式（2）進(jìn)行參數(shù)估計。

3 似然法處理缺失協(xié)變量

3.1 Bayes估計法

用Bayes法進(jìn)行缺失協(xié)變量情形下的參數(shù)估計，基于觀測值（Yi，Xobs，i）估計P(yi|Xi，β) 中的回歸參數(shù)β和P(Xi|α)中的協(xié)變量分布參數(shù)α，設(shè)（α，β）的先驗分布為p(α，β)，損失函數(shù)為平方損失，則參數(shù)的Bayes估計為給定觀測值條件下參數(shù)的后驗期望。參數(shù)的后驗分布函數(shù)：

其 中y=(y1，...，yN) ，Xobs=(Xobs，1，...，Xobs，N) ，Xmis=(Xmis，1，...，Xmis，N)。該后驗分布函數(shù)不容易有顯示表達(dá)式，采用如下步驟獲得參數(shù)的Bayes估計。

第一步：確定參數(shù) (α，β)的聯(lián)合先驗分布p(α，β)，可以假設(shè)p(α，β)=p(α)p(β)即α，β的先驗分布相互獨立，這樣可以對α和β分別假定先驗分布（如定義域上的均勻分布或正態(tài)分布），在Gibbs抽樣時分別從兩個先驗分布中抽取參數(shù)值，也可以假設(shè)p(α，β)≠p(α)p(β)，這樣需要對(α，β)假設(shè)聯(lián)合先驗分布，如二元正態(tài)分布或二維均勻分布，從中抽取參數(shù)值用于后續(xù)Gibbs抽樣。

第三步：Gibbs抽樣以獲取滿足后驗分布p(α，β|y，Xobs)的樣本。抽樣方法如下：

①通過第一步抽取 (α0，β0)；

② 根 據(jù) 式（4）從p(Xmis，i|yi，Xobs，i，α0，β0) 中 抽 取，從而獲得第i個受訪者的完整數(shù)據(jù)集Xmis，i）；

③進(jìn)行第一次迭代，得到參數(shù)β后驗分布p(β|yi，，從中抽取β1；

④獲得參數(shù)α的后驗分布從中抽取α1，這樣就得到第一次迭代后的參數(shù)值(α1,β)；

……

重復(fù)上面的迭代過程，直到（α，β）收斂，這樣就得到一組滿足后驗分布為p(α，β|y，Xobs)的參數(shù)值，記為（α1，β1），重復(fù)m次這一迭代過程，得到m組參數(shù)值（αi，βi），i=1，...，m，參 數(shù) 的 Bayes估 計 量 可 以 用

以上方法假設(shè)協(xié)變量缺失機制為MAR，若缺失機制為NMAR，后驗分布函數(shù)為：

可以容易地將上面的方法推廣到協(xié)變量缺失機制為NMAR的情形。

3.2 EM極大似然估計法

在協(xié)變量含缺失值的情形下，常采用EM算法[5]求參數(shù)的極大似然估計，EM算法是一種迭代算法，每次迭代分兩步進(jìn)行。假設(shè)協(xié)變量的缺失機制為NMAR，首先設(shè)定參數(shù)的迭代初值γ(0)=(α(0)，β(0)，?(0))，設(shè)第 t步迭代后得到的參數(shù)估計值為γ(t)=(α(t)，β(t)，?(t))，第t+1步迭代（t≥0）如下：

E步：根據(jù)（yi，Xi，ri）的聯(lián)合分布p(yi，Xi，ri|α，β，?)=p(yi|Xi，β)p(Xi|α)p(ri|yi，Xi，?)，寫出對數(shù)似然函數(shù)，其中我們最感興趣的參數(shù)是協(xié)變量與因變量之間的回歸系數(shù)β。

由于沒有Xmis，i的信息，需要將其通過積分或求和去除，寫出給定觀測值條件下Xmis，i的密度函數(shù)，求對數(shù)似然 函 數(shù) 關(guān) 于 該 條 件 密 度 的 期 望 ，即|Xobs，i，Yi，ri]，當(dāng)缺失協(xié)變量為連續(xù)型變量時：

若缺失協(xié)變量為離散型，積分換為對缺失協(xié)變量的所有可能值求和，條件分布p(Xmis，i|Xobs，i，yi，ri，γ(t))可以在式（4）基礎(chǔ)上將觀測值ri合并入無缺失觀測值yi,但這時顯式表達(dá)式不容易獲得，嘗試采用Gibbs抽樣的方法獲得缺失 協(xié) 變 量 數(shù) 據(jù) 集 ，具 體 地 ，，采用 Gibbs抽樣的迭代思想，在第 t步參數(shù)γ(t)=(α(t)，β(t)，?(t))已知的條件下，設(shè)Xmis，i為 p1維 向 量，Xobs，i為 p2維 向量，p1+p2=p，從中Gibbs抽樣，其中Xmis，ik為p1維向量Xmis，i的第k個分量。具體步驟如下：

對所有受訪者E得到的結(jié)果為：

一般情況下，我們主要關(guān)心的參數(shù)是反映協(xié)變量與因變量之間的關(guān)系的β，所以α和?為冗余參數(shù)，只需對求最大值得到β(t+1)，它是第t次迭代M步的解，也是第t+1次迭代E步的參數(shù)值。

4 模擬案例

本文通過一個模擬案例比較本文的三種方法（多重插補法、Bayes法、EM極大似然法）與只采用完整數(shù)據(jù)集（刪去含缺失數(shù)據(jù)的受訪者）進(jìn)行回歸分析的完全數(shù)據(jù)法之間的優(yōu)劣。設(shè)樣本量N=250，第i個個體因變量yi為0-1變量，協(xié)變量Xi=(Xi1，Xi2，Xi3),yi通過如下logistic模型獨立生成：

采用如下方式生成協(xié)變量值：

假設(shè)Xi1和Xi2無缺失值，Xi3的缺失機制如下：

可見協(xié)變量的缺失機制為MAR而非MCAR，使用完全數(shù)據(jù)法得到的參數(shù)估計有偏。采用上面的數(shù)據(jù)生成方法，得到1000組樣本量為250的模擬數(shù)據(jù)。在采用上述協(xié)變量數(shù)據(jù)缺失機制，得到1000組含缺失的模擬數(shù)據(jù)。用模擬生成的數(shù)據(jù)對模型（14）進(jìn)行參數(shù)估計，看估計結(jié)果與本文設(shè)定的模型（13）之間的符合程度。

在用Bayes法時，設(shè)先驗分布為無信息先驗π(β)∝1，在多重插補法中，設(shè)插補后得到M=5個完整數(shù)據(jù)集。為獲得基于模擬數(shù)據(jù)參數(shù)估計的標(biāo)準(zhǔn)差（SE）和均方誤差（MSE），對于多重插補法，采用式（2），對于EM極大似然法，采用Ibrahim等（1999）[7]的方法，對于Bayes法，根據(jù)上文所述，通過Gibbs抽樣，在Markov鏈 (αi,βi)收斂后得到m組 (αj，βj),j=1，...，m，取 m=20，可求的標(biāo)準(zhǔn)差和均方誤差。所得結(jié)果如表1所示。

表1 本文三種方法與完全數(shù)據(jù)法參數(shù)估計模擬結(jié)果比較

由以上模擬結(jié)果可見，用完全數(shù)據(jù)法直接刪去協(xié)變量含缺失數(shù)據(jù)受訪者的所有數(shù)據(jù)，僅用完全數(shù)據(jù)樣本得到的結(jié)果與初始模型（13）相比偏差較大，參數(shù)估計的標(biāo)準(zhǔn)差和均方誤差也明顯大于本文其他三種方法所得結(jié)果。所以不推薦使用完全數(shù)據(jù)法進(jìn)行參數(shù)估計。

采用多重插補法、Bayes法和EM極大似然法在協(xié)變量含缺失數(shù)據(jù)情形下進(jìn)行總體參數(shù)估計所得結(jié)果并無明顯優(yōu)劣之分，從無偏性和有效性的兩方面看都可以較好地估計總體參數(shù)，可以根據(jù)實際情況選擇合適的方法進(jìn)行協(xié)變量數(shù)據(jù)缺失情形下的參數(shù)估計。

5 總結(jié)

本文嘗試采用EM算法、Gibbs抽樣法和數(shù)據(jù)擴(kuò)充算法等統(tǒng)計計算方法，采用多重插補法、Bayes法和極大似然法，在缺失機制為MAR或NMAR情形下，估計調(diào)查數(shù)據(jù)中協(xié)變量數(shù)據(jù)缺失情形下，反映協(xié)變量X對因變量y影響大小的回歸系數(shù)。

研究可以進(jìn)一步推廣到協(xié)變量和因變量同時缺失的情形下定義缺失機制、建立似然模型并采用統(tǒng)計算法進(jìn)行回歸系數(shù)估計。還可以推廣到對多個體多時點縱向調(diào)查，在協(xié)變量缺失的情形下，引入統(tǒng)計算法，進(jìn)行回歸系數(shù)估計。