程 豪

(中國科協(xié)創(chuàng)新戰(zhàn)略研究院，北京 100863)

一、引言

眾所周知，線性回歸理論廣泛應用于統(tǒng)計學及其交叉領域。但是，該理論需要滿足較為嚴苛的高斯假設條件，且僅體現(xiàn)平均數(shù)水平，無法全面刻畫因變量在各分位點隨自變量的變化趨勢。相比之下，分位回歸可以處理研究對象異質(zhì)性問題，展示數(shù)據(jù)全貌，有較好的穩(wěn)健性，不要求分布形式。而數(shù)據(jù)的缺失問題普遍存在，如果處理不當，則會得出錯誤的分析結果和研究結論，導致預測與決策的嚴重偏差。

通常，分位回歸中的缺失問題包括因變量缺失和自變量缺失兩類。由于分位回歸不存在似然函數(shù)，目前很多基于似然的處理方法無法直接使用。因此，國內(nèi)外以分位回歸為模型基礎，探索自變量缺失問題仍存在較大的研究空間。

綜上，本文以分位回歸中自變量缺失與因變量有關為選題，旨在通過逆概率加權(Inverse Probability Weighting， IPW)修正現(xiàn)有的多重插補法，提供追溯缺失數(shù)據(jù)、挖掘丟失信息的方法支持，探討中國居民收入影響因素的實際問題。

二、文獻綜述

分位回歸自1978年提出以來，成為一種理論探討和方法應用領域頗具前景的分析工具[1]。數(shù)據(jù)缺失現(xiàn)象普遍存在，一直是備受關注的研究課題。

Hendricks和Koenker利用分位回歸研究了家庭日常用電量與天氣特征之間的關系[2]。Konenker和Machado關注分位數(shù)曲線的數(shù)據(jù)擬合程度評價問題[3]。Bassett和Chen討論了金融市場中多時期的收益問題[4]。Wei等提出了半?yún)?shù)的分位數(shù)模型，并用它建立生長曲線圖[5]。李育安對分位回歸及應用進行簡單介紹[6]。Terry等用分位回歸研究肥胖問題[7]。Wei首次提出了對多元條件分位數(shù)的估計[8]。陳建寶等利用分位回歸對中國居民收入和消費進行了實證分析，但未討論缺失數(shù)據(jù)問題[9]。梅波和田茂再討論了貝葉斯時空分位回歸模型，并用其對北京市PM2.5濃度進行研究[10]。

截止目前，缺失插補方法領域中大多數(shù)情況是完全隨機缺失機制下的完整資料分析法、條件均值插補和基于似然的插補方法[11]。完整資料分析法是一種最簡單的缺失數(shù)據(jù)分析方法。但它功效較低，且導致有偏估計[12]。同樣地，條件均值插補也會帶來有偏估計，因此這些方法不再適用[13]。此時，考慮使用逆概率加權法，在滿足一定分布假定下，基于逆概率加權的估計方程的半?yún)?shù)估計量是有效的、穩(wěn)健的，因此廣受好評[14]。Yi和He對Robins等提出的逆概率加權廣義估計方程進行了拓展，用逆概率加權來修正估計偏差[15-16]。Lipsitz等用逆概率加權方法處理了縱向數(shù)據(jù)缺失問題[17]。但與完整資料分析法類似，逆概率加權法僅僅利用可觀測的完整數(shù)據(jù)信息，造成有效信息的損失。此外，它還會影響估計量的有效性，帶來較大的估計方差。

Wei等提出分位回歸中的多重插補方法，并在此基礎上通過構造壓縮估計量調(diào)整估計偏差[18]。但是，分位回歸中的多重插補法適用于自變量隨機缺失與因變量無關的情況，當自變量隨機缺失與因變量有關時，估計偏差較大。因此，多重插補法有待進一步討論和修正。

三、方法依據(jù)及原理

(一)現(xiàn)有多重插補法的局限

對于分位回歸，缺失插補方法需要解決兩個問題：一是得到x的完整分布，即通過x的條件密度函數(shù)f(x|y，z)生成缺失數(shù)據(jù)的插補值；二是保證或改善估計量的統(tǒng)計性質(zhì)，即在減少估計偏差的同時，有效控制估計方差。

Wei(2012)提出的自變量隨機缺失機制下的多重插補方法，解決了分位回歸的缺失問題，但該方法適用于自變量缺失與因變量無關的情形，當自變量缺失與因變量有關時，該多重插補方法與現(xiàn)有的完整資料分析法一樣，都會帶來較大的估計偏差。

為了解釋多重插補法在自變量缺失與因變量有關時產(chǎn)生估計偏差的原因，令δi是缺失數(shù)據(jù)的示性變量。對于完整資料分析法，需要求解如下估計方程：

∑S[yi，xi，zi，β]*δi=0

(1)

為了得到有效的估計，估計方程的形式為：

E{S[yi，xi，zi，β0]*δi|xi，zi}=0

(2)

其中，β0是參數(shù)真值。

當自變量缺失與因變量無關時，即δi與yi無關時，

E{S[yi，xi，zi，β]*δi}

=E{S[yi，xi，zi，β]|xi，zi}*E{δi}

(3)

因此，E{S[yi，xi，zi，β0]}=0。

但是當自變量缺失與因變量有關時，即δi與yi有關時，

E{S[yi，xi，zi，β]*δi}

≠E{S[yi，xi，zi，β]|xi，zi}*E{δi}

(4)

因此，E{S[yi，xi，zi，β0]}≠0。此時，自變量缺失與因變量有關，多重插補法的參數(shù)估計結果不是無偏估計。

(二)選擇逆概率加權的原因

當自變量缺失與因變量有關時，缺失的數(shù)據(jù)部分與可完整觀測的數(shù)據(jù)部分之間會存在系統(tǒng)差異，如果直接利用可觀測的數(shù)據(jù)信息，依次推斷缺失數(shù)據(jù)，一定會產(chǎn)生推斷結果與真實結果間的差異。換言之，現(xiàn)有的可完整觀測的數(shù)據(jù)點并不能代表總體的分布特征。如果執(zhí)意將現(xiàn)有的每條完整的數(shù)據(jù)記錄以相同的概率納入分析(通常默認為1)，那么可觀測的完整數(shù)據(jù)部分就是總體的一個有偏樣本，估計結果當然也是有偏的。因此，需要對自變量缺失與因變量有關時的每一條可完整觀測的數(shù)據(jù)記錄賦予不同的概率，用以平衡現(xiàn)有可完整觀測數(shù)據(jù)代表總體數(shù)據(jù)分布的偏差。

作為一種減少偏差的修正方法，逆概率加權法最早由Horvitz和Thompson提出，即對每個可觀測的yi的概率取倒數(shù)，作為被觀測的yi的權重，修正由缺失數(shù)據(jù)或有偏抽樣帶來的估計偏差。

(5)

=E{y}>=μ

(6)

可見，利用逆概率加權，總體均值的估計方程一定為無偏估計方程。因此，當自變量缺失與因變量有關時，逆概率加權使得多重插補法的估計方程為無偏估計方程，從而估計出無偏差的估計結果。

(三)用逆概率加權修正多重插補法的原理

作為現(xiàn)有多重插補法的核心內(nèi)容之一，求得x中缺失部分的插補值需要完成概率密度函數(shù)f(x|y，z)的估計，即通過f{y|x，z;β0(τ)}和f{x|z；η(τ)}兩部分完成估計，其中，β0(τ)和η(τ)是兩個真實的分位系數(shù)過程。這是因為，

f(x│y，z)=(f(x，y│z))/(f(y│z))

=f(x|z)f(y|x，z)/f(y│z)

(7)

為了完成逆概率加權對樣本的修正，β0(τ)和η(τ)的估計值一定是逆概率加權以后的參數(shù)估計結果。此外，在完成缺失插補后，逆概率加權多重插補法的分位回歸估計量也是逆概率加權的結果。

(8)

四、模擬研究

(一)基本設定

假定生成模擬數(shù)據(jù)的基礎模型為：

yi=1+xi+zi+(0.5xi+0.5zi)ei

(9)

該模型服從如下假定：

1)自變量向量(xi，zi)服從均值為(4，4)T，方差為(1，1)T的聯(lián)合正態(tài)分布。

2)自變量缺失比例控制在20%左右，缺失概率定義為：P(xi缺失|yi)=1/(1+exp(a+bzi+cyi))。其中，a，b，c為參數(shù)。

3)自變量缺失與因變量有關。通過定義缺失相關參數(shù)(上式中的c)，分別設定c=0.3，0.6，0.9三種情況，表示自變量缺失與因變量有關的不同程度。

4)殘差ei的分布包括兩種假定：標準正態(tài)分布和卡方分布。

此外，模擬過程中包括逆概率加權法(IPW)，多重插補法(MI)，基于觀測概率真值的多重插補法(MI0)，逆概率加權多重插補方法(MI1)四種方法(令MI，MI1，MI0中的m=10)。模擬過程的Monte-Carlo次數(shù)為200，樣本量為1 000。

根據(jù)上述假定，設置如下模擬情形：

情形1：缺失相關參數(shù)為0.3。

情形1-1：缺失相關參數(shù)為0.3，殘差ei服從標準正態(tài)分布。

情形1-2：缺失相關參數(shù)為0.3，殘差ei服從卡方分布。

情形2：缺失相關參數(shù)為0.6。

情形2-1：缺失相關參數(shù)為0.6，殘差ei服從標準正態(tài)分布。

情形2-2：缺失相關參數(shù)為0.6，殘差ei服從卡方分布。

情形3：缺失相關參數(shù)為0.9。

情形3-1：缺失相關參數(shù)為0.9，殘差ei服從標準正態(tài)分布。

情形3-2：缺失相關參數(shù)為0.9，殘差ei服從卡方分布。

(二)模擬結果

1.情形1：缺失相關參數(shù)為0.3

以逆概率加權法(IPW)，多重插補法(MI)，基于觀測概率真值的多重插補法(MI0)三種方法為參照，運行200次Monte-Carlo模擬，實現(xiàn)對逆概率加權多重插補方法(MI1)估計結果的比較研究。表1為缺失相關參數(shù)為0.3時不同分位數(shù)下參數(shù)估計結果(截距相關估計結果已省略)，即對200次Monte-Carlo模擬估計結果取均值、標準差和均方誤差。

表1 缺失相關參數(shù)為0.3時不同分位數(shù)下參數(shù)估計結果(N=1 000)

由表1中關于x和z的估計結果可知，當殘差服從正態(tài)分布時，有如下發(fā)現(xiàn)：1)MI在各分位水平下的估計均值與真值間的差距均為最大，IPW在各分位水平下的估計均值與真值間的差距均為最小。盡管在中分位水平下，MI、MI0和MI1的估計均值與真值間的差距相同，但在低分位水平和高分位水平下，MI0和MI1的估計均值與真值間的差距均小于同等條件下的MI，略大于或等于同等條件下的IPW。2)與其他方法相比，IPW在無偏性上具有明顯的優(yōu)勢。但是，它在有效性和估計精度方面的表現(xiàn)最差。即無論在何種分位水平下，IPW的估計標準差和均方誤差均達到最大值。在高分位水平時表現(xiàn)的尤為明顯。3)MI0和MI1在無偏性、有效性和估計精度上的估計結果均頗為近似。

當殘差服從卡方分布時，關于x和z的估計結果所揭示的規(guī)律有所不同。1)四種方法在低分位水平下的估計均值與真值間的差距相同，IPW在中分位水平下的估計均值與真值間的差距均為最大，其余三種方法在中分位水平下的估計均值與真值間的差距均為相同。在高分位水平下，MI的估計均值與真值間的差距最大，x的系數(shù)估計結果中IPW的估計均值與真值間的差距最小，z的系數(shù)估計結果中MI1的估計均值與真值間的差距最小。2)與殘差服從正態(tài)分布時一致的是，無論在何種分位水平下，IPW的估計標準差和估計精度均達到最大值。在中分位水平和高分位水平時表現(xiàn)的尤為明顯。3)MI0和MI1在無偏性、有效性和估計精度上的估計結果均頗為近似。

2.情形2：缺失相關參數(shù)為0.6

除缺失相關參數(shù)外，其余模擬背景相同。表2是缺失相關參數(shù)為0.6時不同分位數(shù)下參數(shù)估計結果(截距相關估計結果省略)，即對200次Monte-Carlo模擬估計結果取均值、標準差和均方誤差。由表2中關于x和z的估計結果可知，當殘差服從正態(tài)分布時，有如下發(fā)現(xiàn)：1)MI在各分位水平下的估計均值與真值間的差距幾乎均為最大，IPW在各分位水平下的估計均值與真值間的差距均為最小。盡管在低分位水平下，MI，MI0和MI1的估計均值與真值間的差距相同，但在高分位水平下，MI0和MI1的估計均值與真值間的差距均小于同等條件下的MI，略大于或等于同等條件下的IPW。2)與其他方法相比，IPW在無偏性上具有明顯的優(yōu)勢。但是，它在有效性和估計精度方面的表現(xiàn)最差。即無論在何種分位水平下，IPW的估計標準差和均方誤差均達到最大值。3)MI0和MI1在無偏性、有效性和估計精度上的估計結果均頗為近似。

表2 缺失相關參數(shù)為0.6時不同分位數(shù)下參數(shù)估計結果(N=1 000)

當殘差服從卡方分布時，關于x和z的估計結果所揭示的規(guī)律有所不同。1)MI，MI0和MI1三種方法在低分位水平下的估計均值與真值間的差距相同，MI在中分位水平和中分位水平下的估計均值與真值間的差距均為最大。在高分位水平下，四種方法的估計均值與真值間均存在較明顯的差距，x的系數(shù)估計結果中IPW的估計均值與真值間的差距最小，z的系數(shù)估計結果中MI0的估計均值與真值間的差距最小。2)無論在何種分位水平下，IPW的估計標準差和估計精度幾乎均達到最大值(除高分位水平下x系數(shù)的均方誤差達到最大外)。3)MI0和MI1在無偏性、有效性和估計精度上的估計結果均頗為近似。

3.情形2：缺失相關參數(shù)為0.9

表3是缺失相關參數(shù)為0.9時不同分位數(shù)下參數(shù)估計結果(截距相關估計結果已省略)，即對200次Monte-Carlo模擬估計結果取均值、標準差和均方誤差。

由表3中關于x和z的估計結果可知，當殘差服從正態(tài)分布時，有如下發(fā)現(xiàn)：1)MI在各分位水平下的估計均值與真值間的差距幾乎均為最大，IPW在各分位水平下x系數(shù)的估計均值與真值間的差距均為最小。在低分位水平和高分位水平下，MI0和MI1在x系數(shù)的估計均值與真值間的差距均小于同等條件下的MI，略大于或等于同等條件下的IPW。而MI0和MI1在z系數(shù)的估計均值與真值間的差距均最小。2)與其他方法相比，IPW在有效性和估計精度方面的表現(xiàn)最差。除高分位水平下x系數(shù)中MI均方誤差達到最大值以外，IPW的估計標準差和均方誤差均達到最大值。3)MI0和MI1在無偏性、有效性和估計精度上的估計結果均頗為近似。

當殘差服從卡方分布時，關于x和z的估計結果所揭示的規(guī)律有所不同。1)MI、MI0、MI1三種方法在低分位水平下的估計均值與真值間的差距相同，MI在中分位水平下的估計均值與真值間的差距均為最大。在高分位水平下，IPW和MI的估計均值與真值間的差距最大，x的系數(shù)估計結果中IPW的估計均值與真值間的差距最小，z的系數(shù)估計結果中MI1的估計均值與真值間的差距最小。2)除高分位水平下x系數(shù)中MI均方誤差達到最大值以外，IPW的估計標準差和估計精度均達到最大值。3)MI0和MI1在無偏性、有效性和估計精度上的估計結果均頗為近似。

表3 缺失相關參數(shù)為0.9時不同分位數(shù)下參數(shù)估計結果(N=1 000)

(三)可視化研究

為了進一步直觀展示上述模擬結果，進一步研究不同分位水平下四種方法在無偏性和有效性(尤其是無偏性)方面的表現(xiàn)，本文選擇缺失相關系數(shù)為0.3的相關參數(shù)估計結果，繪制圖1。

圖1 參數(shù)估計均值和標準差圖

圖1中，(N，0.1，X)分別表示正態(tài)分布下分位水平為0.1的x系數(shù)估計結果；(N，0.1，Z)分別表示正態(tài)分布下分位水平為0.1的z系數(shù)估計結果；(C，0.1，X)分別表示卡方分布下分位水平為0.1的x系數(shù)估計結果；(C，0.1，Z)分別表示卡方分布下分位水平為0.1的z系數(shù)估計結果，其他同理。圖1突出了不同分位數(shù)下，四種方法無偏性差異的不同規(guī)律。從圖形上來看，隨著分位水平的提高，它們間的差異逐漸變大，高分位水平時最為明顯。為了進一步研究高分位水平下，四種方法的無偏性及變化規(guī)律，本文計算絕對偏差(均值減去真值的絕對值)，并繪制了以不同方法為橫軸、絕對偏差為縱軸的折線圖，見圖2。

圖2 參數(shù)估計偏差圖

圖2中，(N，0.9)分別表示正態(tài)分布下分位水平為0.1的絕對偏差；(C，0.9)分別表示卡方分布下分位水平為0.9的絕對偏差；X.IPW表示IPW對x系數(shù)的估計，其他同理。由圖2可知，無論殘差服從正態(tài)分布還是卡方分布，四種方法按相同順序排列后的絕對偏差整體上呈現(xiàn)相似的變化規(guī)律。相同的是，最大絕對偏差均出現(xiàn)在MI方法處(X.MI和Z.MI)。唯一明顯的不同之處在于最小絕對偏差出現(xiàn)的位置。當殘差服從正態(tài)分布時，x和z的最小絕對偏差均在IPW方法下取得，且優(yōu)勢較為明顯；而當殘差服從卡方分布時，x在IPW、MI0和MI1方法下的最小絕對偏差幾乎處于同一水平，z的最小絕對偏差均在MI1方法下取得，且與MI0頗為接近。

五、應用研究

(一)問題介紹

本文的應用研究基于2010年中國綜合社會調(diào)查(Chinese General Social Survey，CGSS)的部分數(shù)據(jù)，旨在研究年收入的影響因素。研究的對象為1 000名受訪者。本文以年收入(y)為結局變量，周工作時間(x1，含缺失數(shù)據(jù)，缺失率為20%左右)、年齡(x2)、受教育年限(x3)為影響因素?？紤]到這些變量的分布呈現(xiàn)偏態(tài)且均為連續(xù)數(shù)據(jù)，因此可以考慮構建分位回歸模型。模型表達式如下：

Qτ(y)=β0，τ+β1，τx1+β2，τx2+β3，τx3

(10)

(二)研究方法

不難理解，存在缺失數(shù)據(jù)的周工作時間，完整觀測的年齡以及完整觀測的受教育年限都與年收入存在一定的相關關系，即這些自變量在一定程度上會影響或決定了年收入狀況。因此，本文以完整資料分析法(CC)、逆概率加權法(IPW)和多重插補法(MI)為參照，比較逆概率加權多重插補法(MI1)在中國居民收入影響因素問題中的估計結果。其中，多重插補法(MI)和逆概率加權多重插補法(MI1)中m取值為10。

與模擬研究相同，分位回歸中，分位數(shù)水平取值為從0到1、間距為0.02的50個分位數(shù)。整個參數(shù)估計過程包括200次Bootstrap，參數(shù)估計結果包括參數(shù)估計值、200次Bootstrap估計的標準差和參數(shù)顯著性的檢驗P值。

(三)結果解釋

在低分位水平(τ=0.10)，中分位水平(τ=0.50)和高分位水平(τ=0.90)下，不同插補方法的參數(shù)估計結果如表4所示。

表4 不同分位數(shù)水平下的參數(shù)估計結果

由表4可得到一系列方程：

低分位水平：

CC：Qτ(y)=0.27+0.02x1-0.07x2+0.67x3

IPW：Qτ(y)=0.44+0.02x1-0.07x2+0.67x3

MI：Qτ(y)=-0.97+0.02x1-0.06x2+0.71x3

MI1：Qτ(y)=-0.89+0.02x1-0.06x2+0.71x3

中分位水平：

CC：Qτ(y)=4.30+0.01x1-0.06x2+1.14x3

IPW：Qτ(y)=4.11+0.01x1-0.06x2+1.15x3

MI：Qτ(y)=1.85+0.01x1-0.01x2+1.21x3

MI1：Qτ(y)=1.79+0.01x1-0.01x2+1.21x3

高分位水平：

CC：Qτ(y)=-1.68+0.06x1+0.15x2+2.71x3

IPW：Qτ(y)=-1.68+0.06x1+0.15x2+2.71x3

MI：Qτ(y)=-3.54+0.06x1+0.17x2+2.87x3

MI1：Qτ(y)=-4.00+0.07x1+0.17x2+2.87x3

研究表明：分位回歸模型全面展示出，不同年收入水平下，不同人群具有不同的特征和規(guī)律。1)無論在哪種分位水平下，周工作時間(x1)對年收入的貢獻最小、受教育年限(x3)對年收入的貢獻最大，且始終為正向。說明當周工作時間或受教育年限變長時，年度收入會相應增加。其中，受教育年限的影響更為明顯。2)隨著分位水平的升高，受教育年限(x3)對年收入的絕對貢獻(系數(shù)值)有所增加。說明年收入水平越高的人群，具有相對越高的受教育水平。3)年齡隨分位水平的升高，呈現(xiàn)由負變正、不斷增加的趨勢。表明隨著年齡的增加，年收入水平也會相應提高，在一定程度上反映出工齡與收入間存在的關系規(guī)律。綜上所述，受教育年限(或文化程度)成為影響年收入水平的一個關鍵因素。

六、總結與討論

一直以來，提高中國居民工資待遇水平，保障人民生活質(zhì)量都是我們必須面對的重要課題。本文通過模擬研究和應用研究兩部分，將逆概率加權修正后的多重插補法應用于中國居民收入影響因素分析中，嘗試挖掘關鍵因素。

模擬研究表明，在本文涉及的自變量缺失與因變量的相關程度以及殘差服從的分布類型范圍內(nèi)，可以得出如下結論：

第一、逆概率加權法(IPW)對現(xiàn)有多重插補法(MI)的估計偏差進行一定程度的修正，即逆概率加權多重插補法(MI1)比現(xiàn)有多重插補法(MI)在估計偏差上有所改進。

第二、逆概率加權多重插補法(MI1)與基于觀測概率真值的多重插補法(MI0)在無偏性、有效性和估計精度上的估計結果均頗為近似。即在所有分位水平下，逆概率加權多重插補法(MI1)和基于觀測概率真值的多重插補法(MI0)的估計量具有相對較好的無偏性和有效性。

第三、逆概率加權多重插補法(MI1)在不同分位水平下表現(xiàn)出較好的無偏性，有效性和估計精度，可以判定為統(tǒng)計性質(zhì)較佳的估計量。

在模擬中，真實觀測概率是已知的，基于觀測概率真值的多重插補法(MI0)可以使用并作為參照，用于對比研究逆概率加權多重插補法(MI1)。但在實際應用中，真實觀測概率是未知的，基于觀測概率真值的多重插補法(MI0)的討論無法開展，因此，基于觀測概率真值的多重插補法(MI0)僅限于理論層面的探討，逆概率加權多重插補法(MI1)確實對現(xiàn)有多重插補法(MI)有一定程度的改進，但就其程度而言，仍存在著較大的改進空間。將逆概率加權法(IPW)引入現(xiàn)有多重插補法(MI)的思想，有效緩解了當自變量缺失與因變量有關帶來的估計偏差較大的問題，而且并未增加計算復雜度(應用研究中，MI和MI1的單次運行時間分別為11.93秒和11.95秒)。

將逆概率加權多重插補法及其它缺失處理方法應用于中國居民收入影響因素分析中，不難發(fā)現(xiàn)，年度收入會隨著周工作時間、年齡和受教育年限的增長而增加。其中，受教育年限的影響更為明顯。年收入水平較高的人群同時具有相對越高的受教育水平。因此，在中國居民收入的影響因素分析中，受教育程度成為影響收入水平的關鍵因素，相比之下，時間因素對收入水平的影響或貢獻較弱。

綜上所述，逆概率加權多重插補法(MI1)拓寬了現(xiàn)有多重插補法(MI)的適用范圍，確保了估計量具有較好的統(tǒng)計性質(zhì)，而從均方誤差來看，該方法仍然具有可觀的研究前景，使得估計量的性質(zhì)得到更大程度的改善。此外，在中國居民收入影響因素分析中，我們發(fā)現(xiàn)，隨著受教育水平的不斷提高，工作時間因素已不再是制約或影響提高工資待遇的最重要因素，提高學歷水平和文化程度成為改善工資待遇狀況的重要途徑。