協變量數據缺失情形下的參數估計方法

于力超

（中央民族大學 理學院，北京 100081）

0 引言

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

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

本文嘗試引入以上統計計算方法，采用多重插補法、Bayes法和極大似然法，在缺失機制為MAR或NMAR情形下，估計調查數據中協變量數據缺失情形下，協變量關于因變量影響大小的回歸系數。

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

1.1 符號表示

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

1.2 數據缺失機制的概念

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

2 插補法處理缺失協變量

2.1 完全數據法

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

2.2 可用信息法

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

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

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

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

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

2.3 多重插補法

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

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

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

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

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

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

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

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

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

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

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

3 似然法處理缺失協變量

3.1 Bayes估計法

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

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

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

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

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

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

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

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

……

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

以上方法假設協變量缺失機制為MAR，若缺失機制為NMAR，后驗分布函數為：

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

3.2 EM極大似然估計法

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

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

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

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

對所有受訪者E得到的結果為：

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

4 模擬案例

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

采用如下方式生成協變量值：

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

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

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

表1 本文三種方法與完全數據法參數估計模擬結果比較

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

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

5 總結

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

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