樊梅紅， 李婷婷

西南大學 數(shù)學與統(tǒng)計學院，重慶 400715

在統(tǒng)計學中, 有很多方法可以尋找數(shù)據(jù)間的潛在關系, 刻畫數(shù)據(jù)結構． 文獻[1]提出的廣義估計方程(generalized estimating equaiton, GEE)是很常見的一種分析縱向數(shù)據(jù)的統(tǒng)計方法, 在研究數(shù)據(jù)內(nèi)部關系及預測等方面有重要的作用． 此方法的優(yōu)點是即使相關結構被誤判, 所得估計仍然是相合的． 在大數(shù)據(jù)時代, 高維縱向數(shù)據(jù)能比時序數(shù)據(jù)和橫截面數(shù)據(jù)提供更多的信息． 然而數(shù)據(jù)的高維性使模型變得復雜, 降低了模型的估計精度． 帶懲罰項的正則化估計是解決高維數(shù)據(jù)的常用方法． 文獻[2]提出的SCAD懲罰和文獻[3]提出的MCP懲罰是常見的非凸懲罰方法, 具有Oracle性質． 文獻[4-5]將GEE與懲罰函數(shù)相結合, 提出了懲罰廣義估計方程(penalized generalized estimating equation, PGEE), 模擬研究表明該方法在篩選出重要變量的同時得到模型回歸系數(shù)的無偏估計．

在實際應用中, 數(shù)據(jù)往往會呈現(xiàn)異質性． 文獻[6]首次提出分位數(shù)回歸(quantile regression, QR)方法, 可以捕捉整個條件分布的特征． 文獻[7]基于獨立的數(shù)據(jù)結構提出縱向數(shù)據(jù)的線性分位數(shù)回歸模型． 這不可避免地會損失估計效率． 進一步, 文獻[8]考慮縱向數(shù)據(jù)重復觀測樣本間的相關性, 建立分位數(shù)GEE回歸模型, 提高了估計效率． 文獻[9]對縱向數(shù)據(jù)的分位數(shù)回歸模型添加懲罰項, 提出了懲罰分位數(shù)回歸模型．

QR方法對應的損失函數(shù)具有不可微性, 這給數(shù)值計算帶來了很大的難度, 尤其對于高維復雜數(shù)據(jù)來說, 該問題變得更加突出． 受分位數(shù)回歸的啟發(fā), 文獻[10]將分位數(shù)回歸中的非對稱絕對值損失函數(shù)替換為非對稱最小平方損失函數(shù), 提出了期望分位數(shù)(expectile)估計量． Expectile方法不僅繼承了QR方法可以處理異質性的優(yōu)點, 且具有連續(xù)可微的損失函數(shù), 相較QR方法在計算上也有很大的優(yōu)勢． 在獨立同分布的截面數(shù)據(jù)中, 文獻[11-12]將expectile回歸與懲罰函數(shù)相結合, 提出帶有懲罰項的expectile回歸模型, 建立了Oracle性質, 同時實現(xiàn)了變量選擇和異方差識別． 文獻[13]將expectile應用到縱向數(shù)據(jù), 提出了廣義expectile估計方程(generalized expectile estimating equation, GEEE)． 模擬結果顯示, GEEE估計量可以識別出異方差, 在保留分位數(shù)優(yōu)點的同時, 降低了計算難度． 近年來, 作為QR方法的替代, expectile方法受到部分學者的關注, 但在縱向數(shù)據(jù)變量選擇方面的研究還不多見． 本文將截面數(shù)據(jù)的懲罰expectile回歸模型擴展到縱向數(shù)據(jù), 提出PGEEE(penalized generalized expectile estimating equation)估計量． 模擬結果和實證分析顯示, PGEEE估計量不僅可以實現(xiàn)高維數(shù)據(jù)的變量選擇, 并且同時為重要變量的回歸系數(shù)進行估計． 更重要的是, PGEEE方法可以得到一系列τ水平下的變量選擇和模型估計結果, 詳細地刻畫了數(shù)據(jù)的異質結構, 能夠比GEE提供更多的信息．

1 縱向數(shù)據(jù)的expectile回歸和PGEEE模型

1.1 Expectile回歸模型和PGEEE

定義隨機變量Y的τ-expectile值為

(1)

其中τ∈(0, 1),ρτ(θ)=|τ-I(θ≤0)|·θ2是非對稱平方損失函數(shù), I是示性函數(shù)． 由τ-expectile的定義易知, 當τ=0.5時,ρτ(·)等價于經(jīng)典的最小二乘損失函數(shù), 則模型(1)對應經(jīng)典的均值回歸模型,μτ(Y)為隨機變量Y的數(shù)學期望．

假設有縱向樣本數(shù)據(jù)(yij,Xij),i=1,…,n,j=1,…,mi, 滿足如下的expectile線性回歸模型

(2)

yi=Xiβn+εi

(3)

對βn的估計可以通過求解如下目標函數(shù)的最小值來獲得, 即

(4)

考慮重復觀測時個體內(nèi)的相關性, 文獻[13]在縱向數(shù)據(jù)協(xié)變量數(shù)固定的情況下提出了GEEE模型, 即通過求解如下估計方程

(5)

進一步地, 本文在協(xié)變量維數(shù)pn發(fā)散的情況下, 提出縱向數(shù)據(jù)的懲罰非對稱最小二乘PGEEE估計, 即通過求解如下估計方程

Q(βn)=S(βn)-nP′λn(|βn|)Sign(βn)=0

(6)

獲得系數(shù)βn的PGEEE估計． 其中,P′λn(|βn|)=(p′λn(|βn1|), …,p′λn(|βnpn|))T,pλn(t)是一個含有調(diào)節(jié)參數(shù)λn的非負懲罰函數(shù),p′λn(t)為pλn(t)的導數(shù)． Sign(βn)=(sign(βn1), …, sign(βnpn))T, sign(t)=I(t>0)-I(t<0)為符號函數(shù)．P′λn(|βn|)Sign(βn)定義為對應元素相乘得到的向量． 本文考慮MCP和SCAD兩種懲罰方法． MCP懲罰函數(shù)的數(shù)學表達式為

(7)

為簡化模型, 參考文獻[14], 取γ=3． SCAD懲罰函數(shù)的數(shù)學表達式為

(8)

根據(jù)文獻[2]建議取γ=3.7． 此時模型(6)中需要選擇的參數(shù)只有λn, 本文使用BIC準則來選取, 表達式見算法過程．

1.2 求解算法

(10)

Step4: 重復Step2-Step3直至收斂, 并計算λn對應的BIC值, 其表達式為

(11)

其中,df表示λn對應模型所選擇的變量個數(shù)．

2 Oracle性質

注定理1表明所提出的方法可以選出正確的模型, 同時實現(xiàn)對重要變量回歸系數(shù)的參數(shù)估計, 稱為Oracle性質[2]．

定理1的證明:

(12)

成立即可． 根據(jù)表達式, 有

(13)

(βn-βn0)TS(βn)=(βn-βn0)TS(βn0)+(βn-βn0)T[S(βn)-S(βn0)]=I1+I2

(14)

其中

(15)

考慮I11, 有

(16)

(17)

且

(18)

|I1|=Op(pn)‖u‖

(19)

將I2分為兩部分計算, 有

(20)

記

(21)

其中由(A3)知

又因為

(24)

(25)

其中

又

(28)

I2=-Op(pn)‖u‖2

(29)

由(19),(29)式可得, (14)式的值由(29)式控制, 小于0． 易知(13)式中的第二項以nαn2‖u‖ +nbnαn2‖u‖2為界, 因此可以找到一個足夠大的D, 使得(13)式的值完全由(29)式?jīng)Q定． (12)式得證．

(30)

(31)

由(A7)可知, (31)式的符號完全由βj的符號決定． (30)式得證．

即

(32)

(33)

定理證畢．

3 模擬研究

為了研究所提方法的有限樣本性質, 本文比較了不同的懲罰方法及相關結構下所提出方法的效果． 數(shù)據(jù)來源于以下模型

(34)

情形1pn=10,k=9,mi=4,n=50, 100, 200,βn=(-3, 5, 0, 0, 4, 0, 0, 2, 0, 0)T．Ri是參數(shù)為0.9的等相關結構矩陣．

情形2k=2,mi服從參數(shù)為(3, 6)的均勻分布,Ri是參數(shù)為0.9的AR(1)結構矩陣． 其余設置和情形1一樣．

情形3pn=30,n=100, 200．βn=(-3, 5, 0, 0, 4, 0, 0, 2, 0, 0, …, 0)T． 其余設置和情形1一樣．

表1 情形1模擬結果

表2 情形2模擬結果

表3 情形3模擬結果

(i) SCAD和MCP兩種懲罰方法并無明顯的優(yōu)劣之分． FN均為0, 表示所有重要的變量都被識別, FP接近0, 表明噪音變量被選擇的可能性很小;

(ii) 在情形1和情形3中,τ=0.9時, Prob等于1, 而τ=0.5時, Prob的值接近0． 這表明所提出的估計量PGEEE可以在不同的τ水平下, 有效識別出正確的模型, 刻畫數(shù)據(jù)中的異方差結構;

(iii) 在不同的τ水平下, 即使選擇的變量相同, 參數(shù)估計值也可能不同(情形2)． 在此情形下, 估計量的MSE和MAE隨著樣本量增大而減小, 表示該方法可以在識別出異方差的同時實現(xiàn)回歸參數(shù)的一致估計;

(iv) 對比情形1和情形3, 協(xié)變量維數(shù)pn從10增加至30, 結果顯示模型中噪音變量數(shù)量增加時, PGEEE估計表現(xiàn)依然較好, 且估計量MSE減小, 表明該方法可以用于分析高維數(shù)據(jù), 排除無關變量, 識別出重要變量．

(v) 考慮相關結構時估計量的表現(xiàn)總體上優(yōu)于獨立(IND)的情形． 即使相關結構被誤判后, 參數(shù)估計效果依然很好, 尤其使用UN結構時．

4 實證分析

數(shù)據(jù)來自1976年至1982年間對美國經(jīng)濟收入動態(tài)的面板研究, 包含了連續(xù)7年595名民眾的工資水平, 屬于平衡數(shù)據(jù), 更多詳細信息參考文獻[15]． 該研究中, 協(xié)變量包括工作經(jīng)歷E, 工作時間W, 工作職業(yè)O(藍領取1, 否則0), 工作行業(yè)I(制造業(yè)取1, 否則0), 居住地S(居住在南部取1, 否則0), 種族B(黑人取1, 否則0), 是否住在都市統(tǒng)計區(qū)A(如果是取1, 否則0), 是否結婚M(結婚取1, 否則0), 性別F(女性取1, 否則0), 勞動保障U(簽合同取1, 否則0) 及受教育程度D, 響應變量為對數(shù)變換后的工資水平．

表4給出了τ=0.01,0.5,0.95下參數(shù)的PGEEE估計, 其中τ=0． 5對應經(jīng)典的均值回歸估計． 分析結果可知, 不同的懲罰方法和不同的相關結構選出的變量基本一致． 可以看到, 在3個水平下均被選擇的變量有O,B,F,D; 均未被選擇的變量有W． 截距項,B,F的系數(shù)估計隨著τ不同而變化, 圖1a,b為不同種族及性別對應的工資隨時間變化的箱線圖． 男性的工資明顯高于女性, 白人的工資明顯高于黑人． 在τ=0.01時,E被認為是噪音變量, 而在τ=0.5和0.95時被認為是重要變量． 在τ=0.95時, 除了獨立結構下MCP估計外, 工作行業(yè)I, 居住地S, 是否結婚M, 勞動保障U均被剔除在模型外; 而在τ=0.01 和0.5時則被認為是重要變量． 圖1c,d,e,f為這些變量對應的工資分布箱線圖． 以變量S為例, 可以看到, 在低分位點時, 居住在北部的工資要明顯高于南部, 但是在高分位點時, 兩者的區(qū)別并不明顯, 這與PGEEE的估計結果相吻合． 由此可見, 該方法比采用普通最小二乘估計(τ=0.5)挖掘出了更多的信息．

圖1 工資箱線圖

表4 工資數(shù)據(jù)參數(shù)估計結果

5 結語

本文基于expectile提出了高維縱向數(shù)據(jù)的PGEEE估計量, 在實現(xiàn)模型變量選擇的同時, 對模型的回歸系數(shù)進行估計． 在正則條件下本文建立了PGEEE估計量的Oracle性質． 數(shù)值模擬結果顯示, MCP與SCAD懲罰及不同的協(xié)方差結構在變量選擇方面并無明顯差異． 相較于獨立結構, 考慮相關結構時回歸系數(shù)的估計效率更高． 多數(shù)情況下, 不確定結構(UN)的PGEEE估計量具有較好的估計精度． 最后建立工資數(shù)據(jù)的PGEEE模型, 可以看到在不同的τ水平下, 影響工資的因素有所區(qū)別, 同一個因素影響程度也可能不同． 這表明PGEEE可以有效識別數(shù)據(jù)中的異質結構, 比經(jīng)典的懲罰估計方程估計(PGEE)挖掘出更豐富的信息, 更合理地分析了工資的影響因素．