(上海理工大學(xué)理學(xué)院，上海，200093)

1 引言

隨著科學(xué)技術(shù)的不斷發(fā)展，數(shù)據(jù)收集與處理的要求也日益提升.由于數(shù)據(jù)的來(lái)源之廣和數(shù)量的不斷積累，數(shù)據(jù)的結(jié)構(gòu)越來(lái)越復(fù)雜，維數(shù)也越來(lái)越高.如何從海量的數(shù)據(jù)中快速高效的篩選并提取有價(jià)值的信息是我們值得研究和解決的問(wèn)題.在統(tǒng)計(jì)學(xué)中，我們習(xí)慣的稱研究目標(biāo)為響應(yīng)變量，收集到的相關(guān)信息為解釋變量.變量之間的聯(lián)系具有普遍性，一般是借助統(tǒng)計(jì)模型推斷.線性回歸模型是研究變量之間相互關(guān)系的一種有效工具，但當(dāng)統(tǒng)計(jì)模型中有較多的變量時(shí)，如果選擇的變量較少，就會(huì)發(fā)生變量的缺失導(dǎo)致模型會(huì)產(chǎn)生有偏差的估計(jì)和預(yù)測(cè).為了避免這種困境，我們通常會(huì)選擇相對(duì)較多的自變量.但是自變量的過(guò)多又會(huì)產(chǎn)生模型的過(guò)度擬合而影響參數(shù)估計(jì)和預(yù)測(cè)的效率，因此變量選擇是統(tǒng)計(jì)模型和知識(shí)發(fā)現(xiàn)領(lǐng)域非常重要的課題之一.隨著高維數(shù)據(jù)分析需求的日益增加，這一問(wèn)題的研究顯得更加突出而具有現(xiàn)實(shí)意義.

自20世紀(jì)六十年代以來(lái)，變量選擇一直是統(tǒng)計(jì)模型推斷的研究課題之一，許多文獻(xiàn)對(duì)其進(jìn)行了詳細(xì)的研究.傳統(tǒng)的變量選擇有Akaike[1]在熵的概念上提供權(quán)衡估計(jì)模型復(fù)雜度和擬合數(shù)據(jù)優(yōu)良性的標(biāo)準(zhǔn)AIC.從一組可供選擇的模型中選擇最佳模型時(shí)，通常選擇AIC值最小的變量子集.Schwarz[2]提出了BIC的變量選擇方法. BIC準(zhǔn)則與AIC準(zhǔn)則相似，都用于模型選擇，但在懲罰力度上，BIC的懲罰項(xiàng)比AIC的大.在樣本數(shù)量過(guò)多時(shí)，BIC準(zhǔn)則可有效防止模型精度過(guò)高而造成的模型復(fù)雜度過(guò)高.Mallows[3]基于預(yù)測(cè)精度的角度下提出了CP準(zhǔn)則，選擇使CP的值達(dá)到最小的模型作為最優(yōu)模型.Hoerl和Kennard[4]提出了嶺回歸(ridge regression)的變量選擇方法.除了上述方法，常用的傳統(tǒng)變量選擇方法還有逐步回歸、最優(yōu)子集回歸、嶺回歸等.許多研究表明，傳統(tǒng)的變量選擇方法總是由于某些客觀條件的限制，諸如協(xié)變量的維數(shù)很大等而表現(xiàn)的不盡如人意.近年來(lái)基于懲罰函數(shù)類型方法的出現(xiàn)在很大程度上使得變量選擇方法有較大的發(fā)展，開創(chuàng)了變量選擇方法的新景象.Tibshirani[5]提出了LASSO的變量選擇方法.Fan 和Li[6]提出了SCAD的變量選擇方法.Fran 和 Friedman[7]提出了橋回歸(bridge regression)的變量選擇方法.Zhang[8]發(fā)展了MCP的變量選擇方法.Yohai[9]研究了穩(wěn)健的FPE準(zhǔn)則的變量選擇.樊亞莉等[10]作出了縱向數(shù)據(jù)的穩(wěn)健變量選擇方法.Welsh[11]開創(chuàng)了基于bootstrap分層的穩(wěn)健變量選擇.Sophie 等[12]提出了結(jié)合懲罰項(xiàng)的穩(wěn)健變量選擇.

經(jīng)驗(yàn)似然是構(gòu)造置信區(qū)間和檢驗(yàn)非參數(shù)的一門統(tǒng)計(jì)技術(shù)，類似于bootstrap的抽樣特性，這一方法與其它的統(tǒng)計(jì)方法相比有很多突出的優(yōu)點(diǎn).最早，Owen[13]在完全樣本下提出了一種進(jìn)行非參數(shù)類型的統(tǒng)計(jì)推斷方法，并且將經(jīng)驗(yàn)似然方法的應(yīng)用推廣到線性回歸模型的統(tǒng)計(jì)推斷中.Kolaczyk[14]把經(jīng)驗(yàn)似然方法應(yīng)用到了廣義的線性模型.Chen[15-16]研究了非參數(shù)的回歸的經(jīng)驗(yàn)似然統(tǒng)計(jì)推斷.近些年來(lái)，關(guān)于經(jīng)驗(yàn)似然的研究也在不斷的發(fā)展中.Qin 和Lawless[17]研究了約束情況下參數(shù)的估計(jì)方程在經(jīng)驗(yàn)似然中的統(tǒng)計(jì)推斷.還有在復(fù)雜數(shù)據(jù)下所用到的一些經(jīng)驗(yàn)似然方法，如Wang 等[18]在刪失數(shù)據(jù)的一類生存函數(shù)中用到了經(jīng)驗(yàn)似然；Rao等[19]研究了帶缺失數(shù)據(jù)的線性模型的經(jīng)驗(yàn)似然統(tǒng)計(jì)推斷.

本文對(duì)普通的懲罰經(jīng)驗(yàn)似然方法提出了改進(jìn)，詳細(xì)討論基于經(jīng)驗(yàn)似然方法的穩(wěn)健變量選擇問(wèn)題.具體來(lái)說(shuō)，我們把普通經(jīng)驗(yàn)似然約束中的估計(jì)方程穩(wěn)健化，提出了穩(wěn)健經(jīng)驗(yàn)似然估計(jì)，再把懲罰經(jīng)驗(yàn)似然方法和穩(wěn)健估計(jì)方程結(jié)合起來(lái)，提出了一種懲罰穩(wěn)健經(jīng)驗(yàn)似然的壓縮估計(jì)方法，從而達(dá)到穩(wěn)健估計(jì)和變量選擇同時(shí)進(jìn)行.通過(guò)數(shù)值模擬分析，該新方法在變量的準(zhǔn)確度和非零參數(shù)估計(jì)的均方誤差上優(yōu)于普通的懲罰經(jīng)驗(yàn)似然估計(jì).因此，本文的方法在變量選擇等數(shù)據(jù)分析中具有一定的優(yōu)勢(shì).

以下的內(nèi)容安排如下.第2節(jié)介紹本文考慮的模型，以及基于懲罰穩(wěn)健經(jīng)驗(yàn)似然的壓縮估計(jì)方法，第3節(jié)介紹本文的迭代算法和調(diào)節(jié)參數(shù)的選取方法，第4節(jié)進(jìn)行數(shù)值模擬，第5節(jié)是實(shí)證分析，第6節(jié)為本文的總結(jié).

2 模型和方法

2.1 線性回歸模型

考慮簡(jiǎn)單的線性模型：

Y=XTβ+ε,

(2.1)

其中Xi=(Xi1,Xi2,…,Xip)T為解釋變量，i=1,2,…,n，Y=(Y1,Y2,…,Yn)T為響應(yīng)變量，β=(β1,β2,…,βp)T為回歸系數(shù)，ε=(ε1,ε2,…,εn)T為隨機(jī)誤差，ε1,ε2,…,εn獨(dú)立同分布，且E(εi)=0，Var(εi)=σ2，i=1,2,…,n.待估計(jì)的參數(shù)為β，σ2.

2.2 普通的懲罰經(jīng)驗(yàn)似然估計(jì)

普通的最小二乘估計(jì)是基于以下估計(jì)方程得來(lái)的：

(2.2)

根據(jù)Owen[13，20]，模型(2.1)的普通經(jīng)驗(yàn)似然(EL)函數(shù)可表示為

滿足上式的β與pi，可以通過(guò)拉格朗日乘數(shù)法求解.令

(2.3)

普通的經(jīng)驗(yàn)似然函數(shù)為

對(duì)數(shù)經(jīng)驗(yàn)似然函數(shù)為

(2.4)

為了達(dá)到變量選擇的目的，我們?cè)趯?duì)數(shù)經(jīng)驗(yàn)似然函數(shù)的基礎(chǔ)上加懲罰項(xiàng)，得到如下的懲罰經(jīng)驗(yàn)似然函數(shù)

(2.5)

懲罰函數(shù)的選取有多種形式，如文獻(xiàn)[5]提出的LASSO，文獻(xiàn)[8]提出的MCP.本文選取文獻(xiàn)[6]提出的SCAD.SCAD懲罰函數(shù)具有無(wú)偏性、稀疏性、連續(xù)性.無(wú)偏性保證了系數(shù)較大的變量其估計(jì)系數(shù)是漸近無(wú)偏的，這樣避免了不必要的模型偏差.稀疏性保證了將不太重要的系數(shù)壓縮為零，從而起到選擇變量、降低模型復(fù)雜度的效果.連續(xù)性是減少了模型預(yù)測(cè)時(shí)的不穩(wěn)定性，保證了數(shù)據(jù)本身的連續(xù)性.因此本文選取SCAD懲罰函數(shù).

SCAD懲罰函數(shù)具有如下形式：

其中λ>0，a>2.根據(jù)文獻(xiàn)[6]，其一階導(dǎo)數(shù)為

2.3 穩(wěn)健的懲罰經(jīng)驗(yàn)似然估計(jì)

(2.2)式所定義的估計(jì)方法對(duì)數(shù)據(jù)中的異常值是異常敏感的，因此，基于(2.2)式的普通懲罰經(jīng)驗(yàn)似然估計(jì)也會(huì)由于異常值的存在而有較大偏差.為了降低異常值的影響，把(2.2)式穩(wěn)健化，采用有界的得分函數(shù)來(lái)限制估計(jì)方程中異常值的影響.記ψc(x)為對(duì)應(yīng)于Huber函數(shù)的得分函數(shù).函數(shù)ψc(x)的定義為

ψc(x)=min{c,max{-c,x}},

其中調(diào)節(jié)參數(shù)c用來(lái)調(diào)節(jié)估計(jì)的效率和穩(wěn)健性.對(duì)于學(xué)生化殘差Xi,c的取值一般介于1到2之間.

與文獻(xiàn)[21]類似，我們采用如下的權(quán)重函數(shù)來(lái)限制杠桿點(diǎn)的影響：

記

其中σr為σ的絕對(duì)離差估計(jì).

穩(wěn)健經(jīng)驗(yàn)似然函數(shù)可表示為：

Qr(β)的最大值點(diǎn)就是β的穩(wěn)健經(jīng)驗(yàn)似然估計(jì).

加上懲罰函數(shù)的穩(wěn)健經(jīng)驗(yàn)似然函數(shù)為

(2.6)

3 算法

3.1 基于SCAD的穩(wěn)健懲罰經(jīng)驗(yàn)似然估計(jì)的牛頓迭代算法

由于SCAD懲罰函數(shù)是非凸且不可導(dǎo)，直接最小化目標(biāo)函數(shù)并不容易.利用局部二次逼近方法，對(duì)SCAD懲罰函數(shù)在初始值β0展開，有

記

我們使用修改的牛頓迭代算法來(lái)執(zhí)行懲罰經(jīng)驗(yàn)似然函數(shù)的優(yōu)化.具體而言，對(duì)于k=0,1,2,…，我們生成迭代序列

其中

重復(fù)上述迭代過(guò)程，當(dāng)|β(k+1)-β(k)|<ξ時(shí)就停止迭代.這里我們指定ξ=10-6.

3.2 調(diào)節(jié)參數(shù)的選取

為了選擇合適的調(diào)節(jié)參數(shù)，我們對(duì)普通的懲罰經(jīng)驗(yàn)似然估計(jì)和穩(wěn)健的懲罰經(jīng)驗(yàn)似然估計(jì)都運(yùn)用BIC準(zhǔn)則來(lái)選取調(diào)節(jié)參數(shù).為了簡(jiǎn)單起見(jiàn)，在模擬時(shí)令所有的λ相等，即假設(shè)λ(j)=λ,j=1,2,…,p.我們?cè)趨^(qū)間Ω=(0,λmax)中選擇最佳的λ.根據(jù)文獻(xiàn)[23]，

BICλ=log(RSS)+dfλlog(n)/n,

4 數(shù)值模擬

為了比較穩(wěn)健的懲罰經(jīng)驗(yàn)似然方法及壓縮估計(jì)在有限樣本下的性質(zhì),我們通過(guò)數(shù)據(jù)模擬試驗(yàn)來(lái)進(jìn)行說(shuō)明.

4.1 數(shù)據(jù)的產(chǎn)生

我們從模型(2.1)產(chǎn)生模擬數(shù)據(jù)，這里參數(shù)設(shè)定為β=(3,3,0,0,2.5,0,0,0)T，p=8，Xi=(Xi1,Xi2,…,Xip)T的每個(gè)分量服從N(0,1)，樣本大小n=50，對(duì)每組數(shù)據(jù)重復(fù)模擬N=200次.對(duì)于誤差項(xiàng)εi的分布，我們考慮三種情況：

(1)εi服從標(biāo)準(zhǔn)正態(tài)分布N(0,1)，且數(shù)據(jù)無(wú)污染；

(2)εi服從標(biāo)準(zhǔn)正態(tài)分布N(0,1)，但數(shù)據(jù)被污染；

(3)εi服從厚尾的t(2)分布.

4.2 數(shù)據(jù)的3種污染方式

污染1：僅對(duì)解釋變量X的污染, 隨機(jī)污染3%的Xi，并將其替換為Xi-1.2；

污染2：僅對(duì)響應(yīng)變量Y的污染, 隨機(jī)污染3%的Yi，并將其替換為Yi-1.2；

污染3：解釋變量X和響應(yīng)變量Y同時(shí)污染, 隨機(jī)污染3%的Xi，并將其替換為Xi-1.2，隨機(jī)污染3%的Yi，并將其替換為Yi-1.2.

數(shù)值模擬的結(jié)果如下面各表所示.

表1 正態(tài)分布時(shí)數(shù)據(jù)沒(méi)有污染時(shí)的模擬結(jié)果

由表1可知，在模擬數(shù)據(jù)無(wú)污染的情況下，本文所提出的穩(wěn)健懲罰經(jīng)驗(yàn)似然方法及壓縮估計(jì)在非零參數(shù)估計(jì)的變量選擇準(zhǔn)確度上相同于最小二乘估計(jì)和普通的懲罰經(jīng)驗(yàn)似然估計(jì)，但在參數(shù)為零的估計(jì)上明顯優(yōu)于最小二乘估計(jì)，在非零參數(shù)的均方誤差下優(yōu)于普通的懲罰經(jīng)驗(yàn)似然估計(jì).

表2 正態(tài)分布時(shí)數(shù)據(jù)在污染1下的模擬結(jié)果

表3 正態(tài)分布時(shí)數(shù)據(jù)在污染2下的模擬結(jié)果

由表2和表3可知，無(wú)論是在解釋變量污染，還是在響應(yīng)變量污染的情況下，本文所提出的方法在非零參數(shù)估計(jì)的變量選擇準(zhǔn)確度上相同于最小二乘估計(jì)和普通的懲罰經(jīng)驗(yàn)似然估計(jì)，但在參數(shù)為零的估計(jì)上明顯優(yōu)于最小二乘估計(jì)，在非零參數(shù)的均方誤差下優(yōu)于普通的懲罰經(jīng)驗(yàn)似然估計(jì).

表4 正態(tài)分布時(shí)數(shù)據(jù)在污染3下的模擬結(jié)果

由表4可知，即使在解釋變量和相應(yīng)變量都污染的情況下，本文所提出的方法在非零參數(shù)估計(jì)的變量選擇準(zhǔn)確度上相同于最小二乘估計(jì)和普通的懲罰經(jīng)驗(yàn)似然估計(jì)，但在參數(shù)為零的估計(jì)上明顯優(yōu)于最小二乘估計(jì)，在非零參數(shù)的均方誤差下優(yōu)于普通的懲罰經(jīng)驗(yàn)似然估計(jì).

表5 t(2)分布時(shí)數(shù)據(jù)的模擬結(jié)果

由表5可知，不管是模擬數(shù)據(jù)在正態(tài)分布的情況下，還是在自由度為2的厚尾t分布情況下，本文所提出的方法在非零參數(shù)估計(jì)的變量選擇準(zhǔn)確度上仍然相同于最小二乘估計(jì)和普通的懲罰經(jīng)驗(yàn)似然估計(jì)，但在參數(shù)為零的估計(jì)上還是明顯優(yōu)于最小二乘估計(jì)，在非零參數(shù)的均方誤差下優(yōu)于普通的懲罰經(jīng)驗(yàn)似然估計(jì).

在模擬研究中的變量選擇效果上，本文提出的方法相較于普通的懲罰經(jīng)驗(yàn)似然估計(jì)的優(yōu)勢(shì)似乎不太明顯，這可能和參數(shù)β=(3,3,0,0,2.5,0,0,0)T的設(shè)置有關(guān)，因?yàn)槲覀內(nèi)×朔橇愕闹当容^大，信號(hào)比較明顯.

5 實(shí)證分析

在本小節(jié)中，我們將穩(wěn)健的懲罰經(jīng)驗(yàn)似然方法及壓縮估計(jì)應(yīng)用到一組實(shí)際數(shù)據(jù)： Boston房屋數(shù)據(jù).Harrison等[24]用這組數(shù)據(jù)在各種方法下研究了人們對(duì)于干凈空氣的需求，這組數(shù)據(jù)來(lái)源于網(wǎng)址:http://t.cn/RfHTAgY.該數(shù)據(jù)的每個(gè)類的觀察值是均等的，共有506個(gè)觀察值，13個(gè)解釋變量和1個(gè)響應(yīng)變量，其中解釋變量包括：城鎮(zhèn)人犯罪率(CRIM)，住宅用地超過(guò)25000 sq.ft.的比例(ZN)，城鎮(zhèn)非零售商用土地的比例(INDUS)，查理斯河空變量(如果邊界是河流取1，否則取0)(CHAS)，一氧化氮濃度(NOX)，住宅平均房間數(shù)(RM)，1940年之前建成的自用房屋比例(AGE)，到波士頓五個(gè)中心區(qū)域的加權(quán)距離(DIS)，輻射性公路的接近指數(shù)(RAD)，每10000美元的全值財(cái)產(chǎn)稅率(TAX)，城鎮(zhèn)師生比例(PTRATIO),1000(BK-0.63)^2(其中BK指代城鎮(zhèn)中黑人的比例)(B),人口中地位地下者的比例(LSTAT).響應(yīng)變量：自住房的平均房?jī)r(jià)，以千美元計(jì)(MEDV).考慮到實(shí)際的數(shù)據(jù)各項(xiàng)指標(biāo)單位不同的數(shù)值差異較大，因此先將各項(xiàng)指標(biāo)標(biāo)準(zhǔn)化，使其均值為0，方差為1.

為了比較，我們將這組數(shù)據(jù)做了交叉驗(yàn)證來(lái)比較穩(wěn)健的方法和普通的非穩(wěn)健方法的表現(xiàn).每次剔除1個(gè)觀察值，然后用剩下的505個(gè)觀察值的回歸系數(shù)，并用模型在交叉驗(yàn)證過(guò)程中的均方誤差(MSEcv)來(lái)衡量各種方法的優(yōu)劣：

表6 穩(wěn)健與非穩(wěn)健方法對(duì)自住房平均房?jī)r(jià)的MSEcv和Varcv 的估計(jì)比較

MSEcv值越小，表示該變量選擇的效果越好.在表6中可以看出最小二乘估計(jì)的MSEcv最小，但最小二乘估計(jì)不做參數(shù)壓縮和變量選擇，它把所有的參數(shù)值估計(jì)為非零，而穩(wěn)健的方法和非穩(wěn)健的方法都做變量選擇，并且做壓縮估計(jì).綜合來(lái)看，本文提出的穩(wěn)健懲罰方法最好.

6 結(jié)論

本文提出穩(wěn)健的懲罰經(jīng)驗(yàn)似然方法及壓縮估計(jì)，給出了相應(yīng)的算法，并作了大量的模擬研究.根據(jù)數(shù)值模擬研究表明，在數(shù)據(jù)無(wú)污染時(shí)，我們提出的方法在衡量變量選擇的優(yōu)劣中具有較高的準(zhǔn)確度，在非零分量估計(jì)的均方誤差下相對(duì)較??；在數(shù)據(jù)污染的情況下，本文的方法在選擇準(zhǔn)確性和估計(jì)的均方誤差方面的優(yōu)勢(shì)尤為顯著.所以，本文的穩(wěn)健壓縮估計(jì)方法在變量選擇中具有一定可行性的優(yōu)勢(shì).

本論文穩(wěn)健的懲罰經(jīng)驗(yàn)似然方法及壓縮估計(jì)可以推廣到縱向數(shù)據(jù)，也可以結(jié)合其他的懲罰函數(shù)，如LASSO[5]的變量選擇方法，MCP[8]的變量選擇方法等與SCAD的變量選擇方法來(lái)做比較.當(dāng)然，證明本文方法的理論性質(zhì),如變量選擇的相合性，以及壓縮估計(jì)的漸近正態(tài)性，還需要做進(jìn)一步研究.