吳劉倉(cāng)，張家茂，邱貽濤

（昆明理工大學(xué) 理學(xué)院，云南 昆明 650093）

一、引 言

在金融和經(jīng)濟(jì)等領(lǐng)域，經(jīng)常收集到的數(shù)據(jù)不是嚴(yán)格地服從正態(tài)分布，而是服從偏正態(tài)分布，因?yàn)槠珣B(tài)分布能夠較好地刻畫(huà)出數(shù)據(jù)的非對(duì)稱(chēng)性和偏斜程度。一方面，偏態(tài)數(shù)據(jù)是正態(tài)數(shù)據(jù)的進(jìn)一步推廣，是非常常見(jiàn)的一種統(tǒng)計(jì)數(shù)據(jù)；另一方面，在現(xiàn)實(shí)數(shù)據(jù)的采集過(guò)程中，很多抽樣調(diào)查數(shù)據(jù)和實(shí)驗(yàn)數(shù)據(jù)都會(huì)受到無(wú)回答的干擾，或者是因?yàn)槟撤N原因而丟失。因此，對(duì)缺失偏態(tài)數(shù)據(jù)的統(tǒng)計(jì)分析具有重要的理論和實(shí)際意義，筆者旨在研究缺失偏態(tài)數(shù)據(jù)下線(xiàn)性回歸模型參數(shù)的估計(jì)問(wèn)題。

目前，缺失對(duì)稱(chēng)數(shù)據(jù)下回歸模型研究已得到了廣泛關(guān)注，Cheng、Chu等人研究了缺失數(shù)據(jù)下回歸模型中非參數(shù)估計(jì)［1－2］；Wang等人采用回歸插補(bǔ)方法，研究了線(xiàn)性回歸模型中響應(yīng)變量均值的估計(jì)［3］；閆莉等人討論了缺失數(shù)據(jù)下廣義線(xiàn)性模型中參數(shù)置信域問(wèn)題［4］；Little、金勇進(jìn)等人對(duì)缺失數(shù)據(jù)做了詳細(xì)的介紹，并總結(jié)了缺失數(shù)據(jù)的處理方法［5］59－72［6］60－75；Azzalini對(duì)偏正 態(tài)分布 做了大量 研究，并應(yīng)用此模型分析了一組運(yùn)動(dòng)數(shù)據(jù)［7］；Xie等人研究了偏態(tài)數(shù)據(jù)下回歸模型中偏度和尺度參數(shù)的齊次性檢驗(yàn)［8］，而對(duì)缺失偏態(tài)數(shù)據(jù)下回歸模型的研究甚少。

二、偏態(tài)數(shù)據(jù)下線(xiàn)性回歸模型

（一）偏正態(tài)分布

偏正態(tài)分布實(shí)際是一種廣義的正態(tài)分布，一個(gè)隨機(jī)變量Y服從位置參數(shù)μ、尺度參數(shù)σ2和偏度參數(shù)λ的偏正態(tài)分布，記為y～SN（μ，σ2，λ），其密度函數(shù)為：

其中φ（·）、Φ（·）分別是標(biāo)準(zhǔn)正態(tài)分布的概率密度函數(shù)和分布函數(shù)，這種分布的非對(duì)稱(chēng)范圍是（－0．995，0．995）。當(dāng)λ＜0時(shí)，該分布有負(fù)的偏斜；當(dāng)λ＞0時(shí)，該分布有正的偏斜；當(dāng)λ＝0時(shí)，這個(gè)概率密度函數(shù)就是正態(tài)分布的概率密度函數(shù)。

（二）缺失偏態(tài)數(shù)據(jù)下線(xiàn)性回歸模型

偏態(tài)數(shù)據(jù)下線(xiàn)性回歸模型定義如下：

其中xi＝ （xi1，…，xip）T是可以觀(guān)測(cè)的協(xié)變量，β＝（β1，…，βp）Τ是p×1維未知的線(xiàn)性回歸系數(shù)，λ是響應(yīng)變量yi的偏度參數(shù)。設(shè)xi可全部觀(guān)測(cè)到，yi有缺失，δi為指示yi缺失的變量，即：

假定yi滿(mǎn)足隨機(jī)缺失機(jī)制（MAR），即：

其中p（x）表示給定x下y被觀(guān)測(cè)到的概率。

三、完全數(shù)據(jù)下參數(shù)的極大似然估計(jì)

假設(shè)樣本（xi，yi），i＝1，2，…，n來(lái)自模型（2）且獨(dú)立同分布，由式（1）可得對(duì)數(shù)極大似然函數(shù)為：

通常的數(shù)值計(jì)算大多需要使用迭代算法，以下介紹極大似然估計(jì)的常用迭代法。事實(shí)上，這也就是非線(xiàn)性規(guī)劃中求解函數(shù)最大值（最小值）最典型的基本算法，即Gauss－Newton迭代法。

（一）Gauss－Newton迭代法［9］113－114

設(shè)X～f（x，θ），L（θ）＝logf（x，θ），θ∈Θ，則極大似然估計(jì)＝（X）滿(mǎn)足以下必要條件：

在某點(diǎn)θ0處展開(kāi)可得：

因此可視θ0為初值，設(shè)計(jì)以下迭代公式：

其中D（θ）＝［－L″（θ）］－1［L′（θ）］，直到‖θi＋1－θi‖≤ε，ε為預(yù)定的充分小的正數(shù)，如ε＝10－8等，則取θi＋1作為極大似然估計(jì)的近似值。

（二）極大似然估計(jì)的迭代算法

Step1 給定迭代的參數(shù)初值：θ0＝（，，λ0）T。

Step2 給定當(dāng)前值：θi＝（βTi，σ2i，λi）T，迭代θi＋1＝θi＋［－L″（θi）］－1［L′（θi）］。

Step3重復(fù)Step2直到迭代收斂。

四、缺失數(shù)據(jù)下參數(shù)的極大似然估計(jì)

插補(bǔ)方法是處理缺失數(shù)據(jù)的一類(lèi)常用的技術(shù)方法，是指給每一個(gè)缺失數(shù)據(jù)一些替代值，從而得到完整數(shù)據(jù)集；然后使用標(biāo)準(zhǔn)的完全數(shù)據(jù)統(tǒng)計(jì)方法進(jìn)行數(shù)據(jù)分析和推斷。本文在缺失偏態(tài)數(shù)據(jù)下線(xiàn)性回歸模型中，采用不同的插補(bǔ)方法對(duì)缺失數(shù)據(jù)進(jìn)行插補(bǔ)，得到完整數(shù)據(jù)集；再使用完全數(shù)據(jù)下參數(shù)的極大似然估計(jì)對(duì)參數(shù)進(jìn)行估計(jì)。具體插補(bǔ)方法如下：

（一）均值插補(bǔ)（EI）

均值插補(bǔ)法是指對(duì)所有缺失值，用所有觀(guān)測(cè)值的均值進(jìn)行插補(bǔ)。假定可以觀(guān)測(cè)到y(tǒng)1，…，yn1，而yn1＋1，…，yn缺失，即插補(bǔ)值yj為：

得到完全數(shù)據(jù)集，利用完全數(shù)據(jù)下參數(shù)的極大似然估計(jì)，估計(jì)出參數(shù)。

均值插補(bǔ)的特點(diǎn)是操作簡(jiǎn)便，并且對(duì)均值這樣的簡(jiǎn)單變量可以有效地降低其點(diǎn)估計(jì)的偏差，但是由于其缺失值都由均值來(lái)充當(dāng)，因而扭曲了變量的樣本分布，于是均值插補(bǔ)并不適用于偏態(tài)線(xiàn)性回歸模型中回歸系數(shù)、尺度、偏度的參數(shù)估計(jì)。

（二）回歸插補(bǔ)（RI）

回歸插補(bǔ)法是根據(jù)目標(biāo)變量Y和輔助變量X之間的相互關(guān)系建立回歸模型，然后利用已知的輔助變量信息和回歸模型，對(duì)目標(biāo)變量的缺失數(shù)據(jù)進(jìn)行插補(bǔ)的方法。在樣本中，如果變量Y和變量X之間存在高度相關(guān)，可以利用已知數(shù)據(jù)擬合回歸預(yù)測(cè)模型，計(jì)算出回歸替代值，插補(bǔ)出缺失數(shù)據(jù)。

對(duì)于yi～SN（μi，σ2，λ），給定x條件下y的密度函數(shù)為fθ（y｜x），其中θ＝（βT，σ2，λ）T。假定可以觀(guān)測(cè)到y(tǒng)1，…，yn1，而yn1＋1，…，yn缺失，可以利用觀(guān)測(cè)值（x1，y1），…，（xn1，yn1），并采用完全數(shù)據(jù)下參數(shù)的極大似然估計(jì)方法，對(duì)參數(shù)θ進(jìn)行估計(jì)，從而得到。這樣就可以對(duì)缺失值yj（j＝n1＋1，…，n）依照下式進(jìn)行獨(dú)立的參數(shù)隨機(jī)插補(bǔ)：

然后得到插補(bǔ)后的完全數(shù)據(jù)集，利用完全數(shù)據(jù)下參數(shù)的極大似然估計(jì)，得到參數(shù)估計(jì)值。

在偏態(tài)數(shù)據(jù)下線(xiàn)性回歸模型中，由于變量Y和變量X之間存在高度相關(guān)，應(yīng)用回歸插補(bǔ)時(shí)，對(duì)回歸系數(shù)的參數(shù)估計(jì)將會(huì)有十分好的效果，但是對(duì)于相同的xi（i＝1，2，…，n），得到的插補(bǔ)值是相同的，這樣就和均值插補(bǔ)一樣，存在樣本分布扭曲問(wèn)題，將會(huì)造成偏態(tài)數(shù)據(jù)下線(xiàn)性回歸模型中尺度、偏度參數(shù)估計(jì)的不良效果。為了克服這個(gè)缺點(diǎn)，本文在回歸插補(bǔ)方法（RI）的基礎(chǔ)上，針對(duì)缺失偏態(tài)數(shù)據(jù)線(xiàn)性回歸模型（2），提出一種新的迭代插補(bǔ)方法，稱(chēng)之為修正回歸插補(bǔ)（CRI）方法，具體做法如下：

第一步，利用觀(guān)測(cè)值（x1，y1），…，（xn1，yn1），并采用完全數(shù)據(jù)下參數(shù)的極大似然估計(jì)方法，對(duì)參數(shù)θ進(jìn)行估計(jì)，從而得到。

第二步，依次對(duì)第j個(gè)缺失值yj（j＝n1＋1，…，n），依照式（5）進(jìn)行獨(dú)立的參數(shù)隨機(jī)插補(bǔ)。

第三步，將插補(bǔ)值當(dāng)作觀(guān)測(cè)值，得到j(luò)組觀(guān)測(cè)值，重復(fù)第一步，估計(jì)出新的參數(shù)槇θ；重復(fù)第二步，插補(bǔ)出新的缺失值，直到所有缺失值插補(bǔ)完成，即j＝n的時(shí)候 ，將得到參數(shù)的最終估計(jì)。

通過(guò)依次插補(bǔ)缺失值，迭代回歸插補(bǔ)的方法克服了樣本分布扭曲的問(wèn)題，并改善了尺度參數(shù)和偏度參數(shù)的估計(jì)效果。

（三）隨機(jī)回歸插補(bǔ)（RRI）

隨機(jī)回歸插補(bǔ)是在回歸插補(bǔ)的基礎(chǔ)上所作的改進(jìn)，也因考慮到經(jīng)過(guò)回歸后缺失值yj的估計(jì)為yj＝對(duì)于相同的x（i＝1，2，…，n），得到的插補(bǔ)值是i相同的，也存在樣本分布扭曲的問(wèn)題。隨機(jī)回歸插補(bǔ)對(duì)缺失值yj（j＝n1＋1，…，n）依照下式進(jìn)行獨(dú)立的參數(shù)隨機(jī)插補(bǔ)：得到完全數(shù)據(jù)集，利用完全數(shù)據(jù)下參數(shù)的極大似然估計(jì)，估計(jì)出參數(shù)。通過(guò)隨機(jī)地插補(bǔ)缺失值，克服了樣本分布扭曲的缺點(diǎn)，提高了尺度、偏度參數(shù)的估計(jì)效果。

五、Monte Carlo模擬

（一）完全數(shù)據(jù)下的參數(shù)估計(jì)模擬研究

yi，i＝1，2，…，n產(chǎn)生于模型（2），是相互獨(dú)立的隨機(jī)變 量；協(xié) 變 量xi～U（－1，1），取 真 值β＝（－2，3，4）T，σ2＝0．5；在λ＝－0．5、λ＝0、λ＝0．5，而且樣本量為100、200、300時(shí)，模擬1 000次，模擬結(jié)果見(jiàn)表1。

表1 完全數(shù)據(jù)下偏態(tài)線(xiàn)性回歸模型參數(shù)極大似然估計(jì)結(jié)果表

從表1模擬結(jié)果知：隨著樣本量的增加，在不同偏度下，完全偏態(tài)數(shù)據(jù)下線(xiàn)性回歸模型參數(shù)的極大似然估計(jì)的均方誤差（MSE）越來(lái)越小，估計(jì)值越來(lái)越接近真值，說(shuō)明本文的完全偏態(tài)數(shù)據(jù)下線(xiàn)性回歸模型參數(shù)的極大似然估計(jì)效果是良好的，并且不依賴(lài)于偏度，適用于各種偏度情況下的參數(shù)估計(jì)。

（二）缺失數(shù)據(jù)下的參數(shù)估計(jì)模擬研究

yi，i＝1，2，…，n產(chǎn)生于模型（2），是相互獨(dú)立的隨機(jī)變量，協(xié)變量xi～U（－1，1），取真值β＝（2，3，－1）T，σ2＝0．5，λ＝0．5；對(duì)Y分別隨機(jī)缺失5%，10%和30%數(shù)據(jù)后，在插補(bǔ)方法為均值插補(bǔ)（EI）、回歸插補(bǔ)（RI）、隨機(jī)回歸插補(bǔ)（RRI）、修正回歸插補(bǔ)（CRI），而且樣本量n為100，200，300時(shí)，模擬1 000次，其中樣本量為100時(shí)模擬結(jié)果見(jiàn)表2；樣本量為300時(shí)模擬結(jié)果見(jiàn)表3；修正回歸插補(bǔ)方法估計(jì)結(jié)果見(jiàn)表4。

情形1 樣本量n＝100，各種插補(bǔ)方法在不同缺失率下，參數(shù)的估計(jì)結(jié)果比較。

表2 樣本量n＝100時(shí)各種插補(bǔ)方法結(jié)果比較表

情形2 樣本量n＝300，各種插補(bǔ)方法在不同缺失率下，參數(shù)的估計(jì)結(jié)果比較。

表3 樣本量n＝300時(shí)各種插補(bǔ)方法結(jié)果比較表

從表2表3可以得出以下結(jié)論：

1．均值插補(bǔ)后的極大似然參數(shù)估計(jì)隨著缺失率增大，估計(jì)值離真值越來(lái)越遠(yuǎn)，均方誤差（MSE）越來(lái)越大，可見(jiàn)均值插補(bǔ)效果十分差，只適用于缺失率較低情況。

2．回歸插補(bǔ)后的極大似然參數(shù)估計(jì)對(duì)于回歸系數(shù)估計(jì)效果十分好，這與變量Y和輔助變量X之間具有很高的相關(guān)性有關(guān)，但隨著缺失率的增加，尺度參數(shù)估計(jì)值和偏度參數(shù)估計(jì)值離真值越來(lái)越遠(yuǎn)，均方誤差（MSE）逐漸增大，參數(shù)估計(jì)效果比較差。

3．對(duì)比回歸插補(bǔ)，經(jīng)過(guò)隨機(jī)回歸插補(bǔ)后的尺度和偏度參數(shù)估計(jì)，效果有了明顯地改善。

4．經(jīng)過(guò)修正回歸插補(bǔ)后參數(shù)的極大似然參數(shù)估計(jì)，對(duì)回歸系數(shù)、尺度參數(shù)和偏度參數(shù)的估計(jì)效果十分好，而且隨著缺失率的增加，對(duì)所有參數(shù)的估計(jì)都比較穩(wěn)定。參數(shù)估計(jì)效果比隨機(jī)回歸插補(bǔ)后更好，是所有插補(bǔ)方法中參數(shù)估計(jì)總體效果最佳的，十分適合偏態(tài)數(shù)據(jù)下線(xiàn)性回歸模型的參數(shù)估計(jì)。尤其是隨著缺失率和樣本量的增加，以上現(xiàn)象表現(xiàn)得更加明顯，這充分說(shuō)明了筆者提出的修正回歸插補(bǔ)對(duì)缺失偏態(tài)數(shù)據(jù)插補(bǔ)后模型參數(shù)的極大似然估計(jì)，是十分有效的。

情形3 不同樣本量和不同缺失率下，修正回歸插補(bǔ)方法估計(jì)結(jié)果。

從表4可以看出：隨著缺失率的減小，修正回歸插補(bǔ)后的參數(shù)估計(jì)效果越來(lái)越好，符合數(shù)據(jù)缺失下參數(shù)估計(jì)的基本規(guī)律；隨著樣本量的增加，修正回歸插補(bǔ)后的參數(shù)估計(jì)效果越來(lái)越好，進(jìn)一步說(shuō)明了提出的修正回歸插補(bǔ)對(duì)缺失偏態(tài)數(shù)據(jù)插補(bǔ)后模型參數(shù)的極大似然估計(jì)是很好的。

表4 不同樣本量和不同缺失率下修正回歸插補(bǔ)方法估計(jì)結(jié)果表

六、實(shí)例分析

體重指數(shù)（bmi，Y）是與體內(nèi)脂肪總量密切相關(guān)的指標(biāo)，為了簡(jiǎn)單，考慮體重（x1）和性別（x2）兩個(gè)因素，當(dāng)?shù)弥粋€(gè)人的體重和性別就可以簡(jiǎn)要地計(jì)算出體重指數(shù)。該實(shí)例數(shù)據(jù)來(lái)自R軟件sn包中關(guān)于mle的例子，在R軟件中使用sn．mle命令得到估計(jì)參數(shù)，并對(duì)數(shù)據(jù)中心化處理（截距為11．689），處理后的數(shù)據(jù)（bmi）概率密度函數(shù)如圖1。

圖1 bmi概率密度函數(shù)擬合圖

經(jīng)過(guò)處理后的數(shù)據(jù)由圖1可知，體重指數(shù)（bmi）近似服從偏正態(tài)分布，所以令其滿(mǎn)足下列模型：

經(jīng)過(guò)計(jì)算得到完全數(shù)據(jù)下模型參數(shù)估計(jì)如下：

在不同缺失率下對(duì)數(shù)據(jù)隨機(jī)缺失后，利用本文提出的修正回歸插補(bǔ)等方法，計(jì)算得到模型參數(shù)估計(jì)如表5。

從表5可以看出：隨著缺失率的減小，修正回歸插補(bǔ)后的參數(shù)估計(jì)效果越來(lái)越好，本文提出的修正回歸插補(bǔ)方法的表現(xiàn)，是所有插補(bǔ)方法中總體表現(xiàn)最好的。

表5 體重指數(shù)（bmi）的參數(shù)極大似然估計(jì)結(jié)果表

七、結(jié) 論

本文主要目的是研究響應(yīng)變量Y存在偏斜和隨機(jī)缺失下線(xiàn)性回歸模型的參數(shù)估計(jì)問(wèn)題，針對(duì)缺失偏態(tài)數(shù)據(jù)，為了克服樣本分布扭曲缺點(diǎn)，提高模型的回歸系數(shù)、尺度參數(shù)和偏度參數(shù)的估計(jì)效果，提出了一種適合偏態(tài)數(shù)據(jù)下線(xiàn)性回歸模型中缺失數(shù)據(jù)的修正回歸插補(bǔ)方法。通過(guò)隨機(jī)模擬和實(shí)例研究，同均值插補(bǔ)、回歸插補(bǔ)、隨機(jī)回歸插補(bǔ)方法比較，結(jié)果表明所提出的修正回歸插補(bǔ)方法是有用可行的。

［1］ Cheng P E．Nonparametric Estimation of Mean Functionals With Data Missing at andom［J］．J．Amer．Statist Assoc，1994，89（425）．

［2］ Chu C K，Cheng P E．Nonparametric Regression Estimation With Missing Data［J］．Journal of Statist Planning Inference，1995（1）．

［3］ Wang Q H，Rao J N K．Emprical Likelihood for Liner Regression Modles Under Imputation for Missing Responses［J］．Scandinavain Journal of Statistics，2001（4）．

［4］ 閆莉，陳夏．缺失數(shù)據(jù)下廣義線(xiàn)性模型的經(jīng)驗(yàn)似然推斷［J］．統(tǒng)計(jì)與信息論壇，2013（2）．

［5］ Little R J A，Rubin D B．Statistical Analysis With Missing Data［M］．New York：John Wiley ＆ Sons Inc，1987．

［6］ 金勇進(jìn)，邵軍．缺失數(shù)據(jù)的統(tǒng)計(jì)處理［M］．北京：中國(guó)統(tǒng)計(jì)出版社，2009．

［7］ Azzalini A．A Class of Distribution Which Include the Normal Ines［J］．Scandinavain Journal of Statistics，1985（2）．

［8］ Xie F C，Wei B C，Lin J G．Homogeneity Dignostatics for Skew－normal Nonlinear Regression Models［J］．Statistics and Probability Letters，2009（6）．

［9］ 韋博成．參數(shù)統(tǒng)計(jì)教程［M］．北京：高等教育出版社，2006．