丁先文，陳建東，朱小芹

（江蘇理工學(xué)院 數(shù)理學(xué)院，江蘇 常州 213001）

0 引言

關(guān)于線性回歸模型的統(tǒng)計(jì)分析一直都是統(tǒng)計(jì)學(xué)的熱點(diǎn)研究課題?；谧钚《斯烙?jì)方法對回歸模型進(jìn)行統(tǒng)計(jì)分析已被廣泛研究并應(yīng)用在各個(gè)領(lǐng)域。然而，普通的最小二乘回歸（均值回歸）只能描述協(xié)變量對響應(yīng)變量均值的影響，而不能刻畫對響應(yīng)變量條件分布的影響，而且當(dāng)誤差的方差比較大或數(shù)據(jù)中存在異常點(diǎn)時(shí)，最小二乘方法的有效性將備受挑戰(zhàn)。Koenker和Bassett（1978）[1]提出了分位數(shù)回歸模型，隨后對該模型進(jìn)行了系統(tǒng)的研究和推廣。通過估計(jì)不同的條件分位數(shù)函數(shù)，分位數(shù)回歸可以系統(tǒng)地刻畫協(xié)變量對響應(yīng)分布的影響。此外，分位數(shù)回歸模型對誤差分布不作任何假設(shè)，這使得分位數(shù)回歸模型得到了許多研究者的深入研究并在各領(lǐng)域得到了廣泛應(yīng)用。關(guān)于分位數(shù)回歸模型的研究進(jìn)展和詳細(xì)介紹，見Koenker（2005）[2]。

在許多實(shí)際問題中，如抽樣調(diào)查、臨床試驗(yàn)、經(jīng)濟(jì)調(diào)查等，由于被調(diào)查對象不愿回答問題等各種因素，經(jīng)常會(huì)導(dǎo)致缺失數(shù)據(jù)的產(chǎn)生。對缺失數(shù)據(jù)的統(tǒng)計(jì)研究是近年來統(tǒng)計(jì)學(xué)的熱點(diǎn)研究問題。Little和Robin（2014）[3]定義了三種不同的數(shù)據(jù)缺失機(jī)制，即完全隨機(jī)缺失、隨機(jī)缺失（MAR）和不可忽略缺失。在實(shí)際運(yùn)用中，通常假設(shè)數(shù)據(jù)的缺失機(jī)制是MAR，即缺失數(shù)據(jù)只與完全觀測的數(shù)據(jù)有關(guān)。在對缺失數(shù)據(jù)進(jìn)行統(tǒng)計(jì)分析時(shí)，一種常用的方法是只用觀測數(shù)據(jù)進(jìn)行統(tǒng)計(jì)推斷，這在缺失率較大時(shí)會(huì)產(chǎn)生較大的統(tǒng)計(jì)偏差；另一種方法就是對每一個(gè)缺失值采用某種統(tǒng)計(jì)方法進(jìn)行插補(bǔ)，然后對插補(bǔ)后的數(shù)據(jù)集進(jìn)行統(tǒng)計(jì)分析。有關(guān)缺失數(shù)據(jù)方面的介紹，見Little和Robin（2014）[3]。

近年來，對帶有缺失數(shù)據(jù)的分位數(shù)回歸模型的統(tǒng)計(jì)分析引起了一些研究者的興趣。如Sherwood（2013）[4]考慮了協(xié)變量隨機(jī)缺失下分位數(shù)回歸模型的參數(shù)估計(jì)問題；Liu和Yuan（2016）[5]基于經(jīng)驗(yàn)似然方法研究了加權(quán)分位數(shù)回歸的參數(shù)估計(jì)問題。然而，他們的方法都考慮了協(xié)變量隨機(jī)缺失并只利用了觀測數(shù)據(jù)信息進(jìn)行分析，這在數(shù)據(jù)的缺失率較大時(shí)會(huì)帶來較大的估計(jì)誤差。關(guān)于響應(yīng)變量隨機(jī)缺失下分位數(shù)回歸模型的統(tǒng)計(jì)分析很少有文獻(xiàn)報(bào)道。李乃醫(yī)等（2015）[6]研究了響應(yīng)變量隨機(jī)缺失下非線性回歸模型的參數(shù)的經(jīng)驗(yàn)似然置信域問題。本文研究了響應(yīng)變量隨機(jī)缺失下分位數(shù)回歸模型的參數(shù)估計(jì)問題，利用參數(shù)插補(bǔ)方法對缺失的響應(yīng)變量進(jìn)行多重插補(bǔ)，然后基于插補(bǔ)后的數(shù)據(jù)集對回歸模型進(jìn)行參數(shù)估計(jì)。計(jì)算結(jié)果表明該方法在缺失率較大的情況下也可以得到有效的參數(shù)估計(jì)。另外，提出的方法對數(shù)據(jù)中的異常值并不敏感且當(dāng)誤差分布為重尾分布時(shí)，也有較好的估計(jì)結(jié)果。

1 模型

1.1 參數(shù)估計(jì)方法

考慮下面的線性回歸模型：

其中Yi與Xi分別表示響應(yīng)變量及p維協(xié)變量，β是p維的回歸系數(shù)，εi為具有未知分布函數(shù)的隨機(jī)誤差項(xiàng)。在給定Xi的條件下，令Yi的τ條件分位數(shù)為βτ且滿足其中 0＜τ＜1 。設(shè) {Xi,Yi,為來自模型（1）的獨(dú)立同分布的隨機(jī)樣本，其中Xi可被完全觀測，Yi可能缺失。假定δi=1表示Yi可觀測，δi=0表示Yi缺失。本文假定Yi的缺失機(jī)制為隨機(jī)缺失（MAR）,即 πi=P(δi=1|Xi,Yi)=P(δi=1|Xi)，πi在文獻(xiàn)中常被稱為傾向得分（Propensity score）。

對響應(yīng)變量隨機(jī)缺失的情形，參數(shù)估計(jì)可采用成對刪除方法，即：

其中ρτ(t)=t(τ-I(t≤0))為檢查函數(shù)，I(.)為示性函數(shù)。這種估計(jì)方法在缺失率較大的情形下，估計(jì)量的偏差會(huì)變大，使得估計(jì)結(jié)果不可信。處理缺失數(shù)據(jù)的另外一種常見的方法是考慮對每個(gè)缺失的響應(yīng)變量行多次插補(bǔ)，然后基于插補(bǔ)后的數(shù)據(jù)集進(jìn)行參數(shù)估計(jì)。令{y:F(y|Xi)≥u}是Yi在給定協(xié)變量Xi下的第u個(gè)條件分位數(shù)，其中u為來自均勻分布的一個(gè)隨機(jī)數(shù)，F(xiàn)(y|Xi)是給Y在給定協(xié)變量Xi下的條件分布函數(shù)。注意到，在響應(yīng)變量Yi隨機(jī)缺失下，有F(y|X=Xi,δi=1)=F(y|X=Xi,δi=0)成立，這樣就有Qu(yi|X=Xi,δi=1)=Qu(yi|X=Xi,δi=0)。在這里，本文假設(shè)線性模型Qu(Yi|Xi)=XTiβu。在MAR的假設(shè)下，E{δiXi[I(Y＜XTiβu)-u]}=0。

因此，可以得到βu的一個(gè)相合估計(jì)為：

相比于只利用觀測數(shù)據(jù)進(jìn)行參數(shù)估計(jì)，該估計(jì)方法即使在缺失率較大時(shí)也可以得到一個(gè)穩(wěn)健的參數(shù)估計(jì)結(jié)果。另外，與通常的非參數(shù)插補(bǔ)方法相比，即使協(xié)變量維數(shù)p很大時(shí)，該方法也是可行的。

1.2 計(jì)算方法

關(guān)于分位數(shù)回歸模型的參數(shù)估計(jì)，通常的計(jì)算軟件都可實(shí)現(xiàn)。目前較為流行和廣泛采用的方法是利用R軟件中的軟件包quantreg進(jìn)行計(jì)算。假設(shè)為來自模型（1）的獨(dú)立同分布的隨機(jī)樣本，不失一般性，假定前n1個(gè)Yi可觀測，后n-n1個(gè)Yi不可觀測，即將原始樣本分為兩部分。模型的參數(shù)估計(jì)過程如下：

（1）隨機(jī)產(chǎn)生m個(gè)均勻分布的隨機(jī)數(shù){u1,u2,...,um}；

（3）基于觀測數(shù)據(jù){Xk,k=n1+1,...,n}和步驟（2）的結(jié)果，對缺失的Yi進(jìn)行插補(bǔ)，即對每個(gè)Yk，其插補(bǔ)值為

以上計(jì)算過程中的（2）和（4）可以調(diào)用R中的quantreg軟件包實(shí)現(xiàn)。

2 模擬計(jì)算

為實(shí)施模擬，本文從以下模型中產(chǎn)生數(shù)據(jù)：

其中β=(1,2,3)T為三維待估參數(shù)向量，對應(yīng)的Xi的每一個(gè)分量都獨(dú)立同分布于標(biāo)準(zhǔn)正態(tài)分布，Yi根據(jù)模型產(chǎn)生，模型誤差服從以下分布：M1：標(biāo)準(zhǔn)正態(tài)分布N(0,1)；M2：自由度為3的t分布t(3)；M3：混合正態(tài)分布0.1N(0,1)+0.9N(0,10)；M4：混合拉普拉斯分布 0.1Lap(0,1)+0.9Lap(0,10)。假設(shè)缺失概率其中,γ的取值為以下兩種情形：C1：γ=(0.85,1.50,1.50,0.85)T；C2：γ=(-1.80,0.50,1.50,0.85)T，相應(yīng)的響應(yīng)變量Yi的缺失率分別為25%和60%。本文對缺失的Yi進(jìn)行多重插補(bǔ)，插補(bǔ)的次數(shù)設(shè)定為m=10次。通過多次模擬可知，插補(bǔ)10次后，參數(shù)估計(jì)結(jié)果就很穩(wěn)定，并且插補(bǔ)后的估計(jì)結(jié)果對插補(bǔ)次數(shù)并不敏感。將模擬計(jì)算獨(dú)立重復(fù)1000次，計(jì)算結(jié)果如表1和下頁表2所示。表中Bias表示1000次重復(fù)模擬的參數(shù)估計(jì)的均值與真實(shí)值的絕對偏差，SD表示1000次模擬的參數(shù)估計(jì)的標(biāo)準(zhǔn)差，RMS表示1000次模擬的參數(shù)估計(jì)的均值與真實(shí)值之差的平方和的平方根。

表1 缺失率為25%時(shí)的模擬計(jì)算結(jié)果

表2 缺失率為60%時(shí)的模擬計(jì)算結(jié)果

從表1和表2可以看出：

（1）對兩種估計(jì)方法，減少缺失率有利于減小估計(jì)的偏差和提高估計(jì)的精度；

（2）針對不同的誤差M1至M4，基于文中提出的插補(bǔ)方法得到的參數(shù)估計(jì)結(jié)果比基于完全觀測數(shù)據(jù)得到的結(jié)果具有較小的Bias、SD和RMS，這說明提出的多重插補(bǔ)方法可以減小估計(jì)偏差并能提高估計(jì)的有效性；

（3）基于插補(bǔ)方法得到的參數(shù)估計(jì)結(jié)果在不同的分位點(diǎn)處表現(xiàn)都很好。

3 總結(jié)

本文研究了在響應(yīng)變量隨機(jī)缺失下分位數(shù)回歸模型的參數(shù)估計(jì)問題。傳統(tǒng)的基于成對刪除數(shù)據(jù)的估計(jì)方法沒有利用所有可觀測到的協(xié)變量的信息，在數(shù)據(jù)的缺失率較大時(shí)容易產(chǎn)生較大的偏差，降低了估計(jì)效率。本文首先基于觀測數(shù)據(jù)得到了模型的參數(shù)估計(jì)；其次在MAR假定下，對缺失的響應(yīng)變量進(jìn)行了多重插補(bǔ)，從而得到了插補(bǔ)后的數(shù)據(jù)集；最后基于插補(bǔ)后的數(shù)據(jù)集對分位數(shù)回歸模型進(jìn)行參數(shù)估計(jì)。該插補(bǔ)方法在數(shù)據(jù)的缺失率較大時(shí)依然有效，并且由于采用的是參數(shù)插補(bǔ)方法，即使當(dāng)協(xié)變量維數(shù)很高時(shí)，方法依然有效。模擬計(jì)算說明了該方法的有效性。本文的方法可以應(yīng)用于微觀經(jīng)濟(jì)、醫(yī)藥追蹤試驗(yàn)和抽樣調(diào)查等帶有缺失數(shù)據(jù)的各種領(lǐng)域。

參考文獻(xiàn)：

[1]Koenker R,Bassett Jr G.Regression Quantiles[J].Econometrica：Journal of the Econometric Society,1978.

[2]Koenker R.Quantile Regression[M].Cambridge：Cambridge University Press,2005.

[3]Little R J A,Rubin D B.Statistical Analysis with Missing Data[M].New Jersey：John Wiley&Sons,2014.

[4]Sherwood B,Wang L,Zhou X H.Weighted Quantile Regression for Analyzing Health Care Cost Data With Missing Covariates[J].Statistics in Medicine,2013,32(28).

[5]Liu T,Yuan X.Weighted Quantile Regression With Missing Covariates Using Empirical likelihood[J].Statistics,2016,50(1).

[6]李乃醫(yī),李永明,韋盛學(xué).缺失數(shù)據(jù)下非線性分位數(shù)回歸模型的光滑經(jīng)驗(yàn)似然推斷[J].統(tǒng)計(jì)與決策,2015，(1).