變系數(shù)模型的穩(wěn)健變量選擇與結(jié)構(gòu)識(shí)別

聶黎雯

(重慶師范大學(xué) 數(shù)學(xué)科學(xué)學(xué)院，重慶 401331)

0 引言

變系數(shù)模型(Varying Coefficient Models，簡(jiǎn)稱“VCM”)是Hastie和Tibshirani在1993年提出來(lái)的一類半?yún)?shù)模型，一般形式為

Y=XTβ(U)+ε，

(1)

其中：Y∈R為模型中的響應(yīng)變量；X=(X1,…,Xp)T∈RP是p維解釋變量；U為指標(biāo)變量；β(·)=(β1(·),…,βp(·))T是模型中p×1維的函數(shù)向量；ε為模型誤差.

目前，眾多學(xué)者對(duì)該模型的變量選擇和結(jié)構(gòu)識(shí)別進(jìn)行了研究.Lv[1]基于雙重SCAD懲罰函數(shù)和秩回歸，提出了能同時(shí)解決參數(shù)估計(jì)和結(jié)構(gòu)識(shí)別的方法.Ma等[2]基于核密度估計(jì)與SCAD懲罰函數(shù)解決了該模型的統(tǒng)一變量選擇問(wèn)題.Tang等[3]基于自適應(yīng)LASSO懲罰和兩步迭代程序?qū)崿F(xiàn)了變量選擇和結(jié)構(gòu)識(shí)別.Zhao[4]基于分位數(shù)估計(jì)以及核函數(shù)實(shí)現(xiàn)了變系數(shù)模型的統(tǒng)一變量選擇.雖然后面兩篇文章能同時(shí)解決變量選擇和結(jié)構(gòu)識(shí)別，但是基于最小二乘法和分位數(shù)回歸.最小二乘法不是穩(wěn)健估計(jì)，分位數(shù)回歸會(huì)因分位數(shù)點(diǎn)不同導(dǎo)致結(jié)果不同.Zou等[5]首次提出了綜合各個(gè)分位數(shù)點(diǎn)的復(fù)合分位數(shù)回歸(CQR).在異方差的情況下，Yang等[6]考慮了該模型基于復(fù)合分位數(shù)回歸的穩(wěn)健變量選擇.Guo等[7]基于加權(quán)復(fù)合分位數(shù)估計(jì)與SCAD懲罰函數(shù)實(shí)現(xiàn)了該模型的變量選擇.Matthew Pietrosanu等[8]討論了變系數(shù)模型的分位數(shù)回歸和復(fù)合分位數(shù)回歸與組Lasso相結(jié)合的變量選擇問(wèn)題.但上述研究沒(méi)有考慮常系數(shù)與變系數(shù)的區(qū)分.因此，本文基于復(fù)合分位數(shù)的穩(wěn)健性與有效性，結(jié)合SCAD懲罰函數(shù)和B樣條基函數(shù)近似提出了一種新的雙懲罰穩(wěn)健估計(jì)程序，它能夠同時(shí)進(jìn)行變量選擇和結(jié)構(gòu)識(shí)別，而且所得估計(jì)具有更好的穩(wěn)健性和有效性.

1 模型的建立

假設(shè)(Yi,Xi,Ui)，i=1,…,n為來(lái)自變系數(shù)模型(1)的一組獨(dú)立同分布樣本.利用B樣條基函數(shù)去近似模型(1)中的β(·).令B(·)=(B1(·),…,Bqn(·))T為m階的B樣條基，則每個(gè)函數(shù)系數(shù)βl(·)都可以用如下表達(dá)式近似：

(2)

其中：γl=(γl,1,…,γl,qn)T是系數(shù)向量；qn=kn+m+1是基函數(shù)個(gè)數(shù)；kn為內(nèi)節(jié)點(diǎn)的個(gè)數(shù).基于以上近似，則模型(1)可寫為

(3)

令πi=Ip?B(Ui)·Xi，Ip是p×p的單位矩陣.模型(3)可被改寫為

(4)

則變系數(shù)模型的復(fù)合分位數(shù)函數(shù)為

(5)

ρτk(x)=τkxI(x≥0)-(1-τk)xI(x<0)，

(6)

其中：τk=k/(q+1)；q為所取分位數(shù)點(diǎn)的個(gè)數(shù)；cτk是隨機(jī)誤差ε的τk分位點(diǎn).

本文的目的是從模型(5)中識(shí)別出常系數(shù)變量和變系數(shù)變量以及系數(shù)為零的變量.設(shè)pλ′(·)為SCAD懲罰函數(shù)(Fan和Li[9])的一階導(dǎo)函數(shù)：

(7)

其中，a>2，θ>0并且pλ(0)=0以及非負(fù)的懲罰參數(shù)λ.為了進(jìn)行變量選擇和參數(shù)識(shí)別，將以下目標(biāo)函數(shù)最小化，

(8)

其中：

對(duì)于目標(biāo)函數(shù)(8)來(lái)說(shuō)，第一個(gè)懲罰pλ1(·)是用來(lái)控制模型的稀疏程度.顯然，如果‖γl‖2被壓縮為0，則βl(·)=0.從第二個(gè)懲罰函數(shù)可以看出

其中

fl={βl(U2)-βl(U1),…,
βl(Un)-βl(Un-1)}T.

2 算法與調(diào)節(jié)參數(shù)的選取

2.1 算法

首先，本文采用MM算法，構(gòu)造一個(gè)優(yōu)化函數(shù)去逼近目標(biāo)函數(shù)L(γ)，然后最小化優(yōu)化函數(shù).以L(γ)為基礎(chǔ)，函數(shù)(5)能夠被優(yōu)化為

(9)

進(jìn)一步，除去常數(shù)部分的目標(biāo)函數(shù)Ln(γ)可以被近似替代為

(10)

為了提高局部二次近似方法的精確性，將

替換為

以及

替換為

其中:?為擾動(dòng)項(xiàng)(?>0).關(guān)于?取值可以參考文獻(xiàn)[10].進(jìn)一步利用牛頓迭代算法可以獲得(10)的數(shù)值解.

2.2 參數(shù)的選取

為了實(shí)現(xiàn)這個(gè)算法，需要選擇B樣條的階數(shù)m，內(nèi)結(jié)點(diǎn)的個(gè)數(shù)kn以及兩個(gè)懲罰參數(shù)λ1，λ2.對(duì)于變系數(shù)模型，通常使用低階樣條，所以本文使用三階B樣條(m=3)，由于B樣條內(nèi)結(jié)點(diǎn)的個(gè)數(shù)對(duì)參數(shù)估計(jì)不敏感，因此，本文取kn=[n1/(2m+3)]，其中[a]表示取值不超過(guò)a的最大整數(shù).通過(guò)最小化BIC對(duì)λ=(λ1,λ2)進(jìn)行選擇，BIC的定義為

(11)

3 數(shù)據(jù)模擬

例1數(shù)據(jù)來(lái)源于以下部分線性變系數(shù)模型：

考慮以下3種誤差分布：(1)標(biāo)準(zhǔn)正態(tài)分布N(0,1);(2)自由度為3的t分布t3;(3)Lognormal (0,1)分布.在模擬中，考慮n=200,300，重復(fù)模擬200次.由于復(fù)合分位數(shù)估計(jì)不依賴q的選取，所以本文考慮q=9.同時(shí)在數(shù)值模擬中對(duì)比以下兩種方法：(1)懲罰最小二乘估計(jì)(PLS)；(2)懲罰分位數(shù)回歸(PQR).考慮分位數(shù)τ=0.5，0.75，相應(yīng)的估計(jì)記為PQR0.5、PQR0.75.同時(shí)為了說(shuō)明PCQR估計(jì)的有效性，采用以下幾個(gè)評(píng)價(jià)指標(biāo)：

(1)均方誤差平方根(Square root of average square errors,RASE)，其定義為

(2)正確識(shí)別零系數(shù)的平均個(gè)數(shù)：C.

(3)非零系數(shù)識(shí)別為零系數(shù)的平均個(gè)數(shù)：IC.

(4)正確識(shí)別變系數(shù)的平均個(gè)數(shù)：CV.

(5)正確識(shí)別常系數(shù)的平均個(gè)數(shù)：CC.

(6)正確識(shí)別真實(shí)模型的比率：CF.

(7)真實(shí)模型下估計(jì)的RASE：Oracle.

(8)懲罰估計(jì)的RASE：Penalized.

具體的數(shù)值模擬結(jié)果如表1所列.

從表1中可以得出以下結(jié)論：首先，當(dāng)n=200，誤差服從正態(tài)分布時(shí)，PLS的Oracle估計(jì)雖然具有最小的RASE，但是PCQR估計(jì)的CF是3種方法中表現(xiàn)最好的.當(dāng)誤差分布服從厚尾分布(t3、Lognormal)時(shí)，發(fā)現(xiàn)盡管PQR估計(jì)整體上都比PLS估計(jì)好，但是關(guān)于CF和RASE方面，PCQR估計(jì)表現(xiàn)最好，這也就意味著PCQR估計(jì)不僅比PLS估計(jì)穩(wěn)健同時(shí)還比PQR估計(jì)有效.最后，隨著變量的增加，PCQR估計(jì)選出真實(shí)模型的比率明顯高于PLS和PQR且越來(lái)越來(lái)接近1.同時(shí)PCQR的Oracle估計(jì)具有最小的RASE，且懲罰下的RASE與 Oracle估計(jì)的RASE差距最小，這也表明本文的PCQR具有Oracle性質(zhì).

表1 例1的數(shù)值結(jié)果

4 實(shí)例分析

為了更好地驗(yàn)證PCQR方法的有效性，本節(jié)將所提方法運(yùn)用到Boston住房數(shù)據(jù)中.根據(jù)Wang和Xia(2009)[11]的建議，將房子價(jià)格的中位數(shù)(medv)作為響應(yīng)變量Y，地位較低人口占總?cè)丝诘陌俜致?lstat)作為指標(biāo)變量U，再將以下6個(gè)變量作為協(xié)變量X：鎮(zhèn)上人均犯罪率(crim)、氮氧化合物濃度(nox)、每座住宅的房間數(shù)的均值(rm)以及在1940年之前建造的房屋所占比率(age)、房屋全值物業(yè)稅率(tax)、鎮(zhèn)上學(xué)生人數(shù)和老師人數(shù)的比率(ptratio).先將所有的協(xié)變量都看作變系數(shù)，用本文所提出的方法進(jìn)行選擇和識(shí)別.同時(shí)，將PCQR與PLS和PQR進(jìn)行對(duì)比.為了檢驗(yàn)?zāi)Ｐ偷念A(yù)測(cè)能力，將506個(gè)數(shù)據(jù)隨機(jī)分為訓(xùn)練集(training data)和預(yù)測(cè)集(testing data)，其中300個(gè)數(shù)據(jù)作為訓(xùn)練集用于變量選擇和參數(shù)估計(jì)，剩余的206個(gè)數(shù)據(jù)則看成是預(yù)測(cè)集用來(lái)評(píng)估模型的預(yù)測(cè)能力.本文利用MAPE(mean absolute prediction error)來(lái)衡量模型的外預(yù)測(cè)能力，其定義為

實(shí)驗(yàn)結(jié)果如表2所列.

重復(fù)模擬200次，用從6個(gè)變量中做變量選擇后的平均個(gè)數(shù)(Size)刻畫模型的復(fù)雜程度.若擬合效果越好，預(yù)測(cè)精度越高，則外預(yù)測(cè)(MAPE)越小.從表2的結(jié)果可以看出，PCQR的MAPE明顯小于其余兩種估計(jì)方法的MAPE，且模型的復(fù)雜程度低于其余兩種方法，說(shuō)明本文提出的估計(jì)方法的有效性.

表2 基于200次重復(fù)，各種方法的MAPE

5 結(jié)語(yǔ)

本文針對(duì)變系數(shù)模型變量選擇和結(jié)構(gòu)識(shí)別問(wèn)題，提出了能夠同時(shí)對(duì)模型進(jìn)行變量選擇和結(jié)構(gòu)識(shí)別的新方法.通過(guò)數(shù)值模擬分析，PCQR比現(xiàn)有方法更加穩(wěn)健和有效.同時(shí)實(shí)證分析也說(shuō)明PCQR比現(xiàn)有方法有更好的預(yù)測(cè)精度.