樊梅紅, 李婷婷
西南大學(xué) 數(shù)學(xué)與統(tǒng)計(jì)學(xué)院,重慶 400715
在統(tǒng)計(jì)學(xué)中, 有很多方法可以尋找數(shù)據(jù)間的潛在關(guān)系, 刻畫數(shù)據(jù)結(jié)構(gòu). 文獻(xiàn)[1]提出的廣義估計(jì)方程(generalized estimating equaiton, GEE)是很常見的一種分析縱向數(shù)據(jù)的統(tǒng)計(jì)方法, 在研究數(shù)據(jù)內(nèi)部關(guān)系及預(yù)測等方面有重要的作用. 此方法的優(yōu)點(diǎn)是即使相關(guān)結(jié)構(gòu)被誤判, 所得估計(jì)仍然是相合的. 在大數(shù)據(jù)時代, 高維縱向數(shù)據(jù)能比時序數(shù)據(jù)和橫截面數(shù)據(jù)提供更多的信息. 然而數(shù)據(jù)的高維性使模型變得復(fù)雜, 降低了模型的估計(jì)精度. 帶懲罰項(xiàng)的正則化估計(jì)是解決高維數(shù)據(jù)的常用方法. 文獻(xiàn)[2]提出的SCAD懲罰和文獻(xiàn)[3]提出的MCP懲罰是常見的非凸懲罰方法, 具有Oracle性質(zhì). 文獻(xiàn)[4-5]將GEE與懲罰函數(shù)相結(jié)合, 提出了懲罰廣義估計(jì)方程(penalized generalized estimating equation, PGEE), 模擬研究表明該方法在篩選出重要變量的同時得到模型回歸系數(shù)的無偏估計(jì).
在實(shí)際應(yīng)用中, 數(shù)據(jù)往往會呈現(xiàn)異質(zhì)性. 文獻(xiàn)[6]首次提出分位數(shù)回歸(quantile regression, QR)方法, 可以捕捉整個條件分布的特征. 文獻(xiàn)[7]基于獨(dú)立的數(shù)據(jù)結(jié)構(gòu)提出縱向數(shù)據(jù)的線性分位數(shù)回歸模型. 這不可避免地會損失估計(jì)效率. 進(jìn)一步, 文獻(xiàn)[8]考慮縱向數(shù)據(jù)重復(fù)觀測樣本間的相關(guān)性, 建立分位數(shù)GEE回歸模型, 提高了估計(jì)效率. 文獻(xiàn)[9]對縱向數(shù)據(jù)的分位數(shù)回歸模型添加懲罰項(xiàng), 提出了懲罰分位數(shù)回歸模型.
QR方法對應(yīng)的損失函數(shù)具有不可微性, 這給數(shù)值計(jì)算帶來了很大的難度, 尤其對于高維復(fù)雜數(shù)據(jù)來說, 該問題變得更加突出. 受分位數(shù)回歸的啟發(fā), 文獻(xiàn)[10]將分位數(shù)回歸中的非對稱絕對值損失函數(shù)替換為非對稱最小平方損失函數(shù), 提出了期望分位數(shù)(expectile)估計(jì)量. Expectile方法不僅繼承了QR方法可以處理異質(zhì)性的優(yōu)點(diǎn), 且具有連續(xù)可微的損失函數(shù), 相較QR方法在計(jì)算上也有很大的優(yōu)勢. 在獨(dú)立同分布的截面數(shù)據(jù)中, 文獻(xiàn)[11-12]將expectile回歸與懲罰函數(shù)相結(jié)合, 提出帶有懲罰項(xiàng)的expectile回歸模型, 建立了Oracle性質(zhì), 同時實(shí)現(xiàn)了變量選擇和異方差識別. 文獻(xiàn)[13]將expectile應(yīng)用到縱向數(shù)據(jù), 提出了廣義expectile估計(jì)方程(generalized expectile estimating equation, GEEE). 模擬結(jié)果顯示, GEEE估計(jì)量可以識別出異方差, 在保留分位數(shù)優(yōu)點(diǎn)的同時, 降低了計(jì)算難度. 近年來, 作為QR方法的替代, expectile方法受到部分學(xué)者的關(guān)注, 但在縱向數(shù)據(jù)變量選擇方面的研究還不多見. 本文將截面數(shù)據(jù)的懲罰expectile回歸模型擴(kuò)展到縱向數(shù)據(jù), 提出PGEEE(penalized generalized expectile estimating equation)估計(jì)量. 模擬結(jié)果和實(shí)證分析顯示, PGEEE估計(jì)量不僅可以實(shí)現(xiàn)高維數(shù)據(jù)的變量選擇, 并且同時為重要變量的回歸系數(shù)進(jìn)行估計(jì). 更重要的是, PGEEE方法可以得到一系列τ水平下的變量選擇和模型估計(jì)結(jié)果, 詳細(xì)地刻畫了數(shù)據(jù)的異質(zhì)結(jié)構(gòu), 能夠比GEE提供更多的信息.
定義隨機(jī)變量Y的τ-expectile值為
(1)
其中τ∈(0, 1),ρτ(θ)=|τ-I(θ≤0)|·θ2是非對稱平方損失函數(shù), I是示性函數(shù). 由τ-expectile的定義易知, 當(dāng)τ=0.5時,ρτ(·)等價于經(jīng)典的最小二乘損失函數(shù), 則模型(1)對應(yīng)經(jīng)典的均值回歸模型,μτ(Y)為隨機(jī)變量Y的數(shù)學(xué)期望.
假設(shè)有縱向樣本數(shù)據(jù)(yij,Xij),i=1,…,n,j=1,…,mi, 滿足如下的expectile線性回歸模型
(2)
yi=Xiβn+εi
(3)
對βn的估計(jì)可以通過求解如下目標(biāo)函數(shù)的最小值來獲得, 即
(4)
考慮重復(fù)觀測時個體內(nèi)的相關(guān)性, 文獻(xiàn)[13]在縱向數(shù)據(jù)協(xié)變量數(shù)固定的情況下提出了GEEE模型, 即通過求解如下估計(jì)方程
(5)
進(jìn)一步地, 本文在協(xié)變量維數(shù)pn發(fā)散的情況下, 提出縱向數(shù)據(jù)的懲罰非對稱最小二乘PGEEE估計(jì), 即通過求解如下估計(jì)方程
Q(βn)=S(βn)-nP′λn(|βn|)Sign(βn)=0
(6)
獲得系數(shù)βn的PGEEE估計(jì). 其中,P′λn(|βn|)=(p′λn(|βn1|), …,p′λn(|βnpn|))T,pλn(t)是一個含有調(diào)節(jié)參數(shù)λn的非負(fù)懲罰函數(shù),p′λn(t)為pλn(t)的導(dǎo)數(shù). Sign(βn)=(sign(βn1), …, sign(βnpn))T, sign(t)=I(t>0)-I(t<0)為符號函數(shù).P′λn(|βn|)Sign(βn)定義為對應(yīng)元素相乘得到的向量. 本文考慮MCP和SCAD兩種懲罰方法. MCP懲罰函數(shù)的數(shù)學(xué)表達(dá)式為
(7)
為簡化模型, 參考文獻(xiàn)[14], 取γ=3. SCAD懲罰函數(shù)的數(shù)學(xué)表達(dá)式為
(8)
根據(jù)文獻(xiàn)[2]建議取γ=3.7. 此時模型(6)中需要選擇的參數(shù)只有λn, 本文使用BIC準(zhǔn)則來選取, 表達(dá)式見算法過程.
(10)
Step4: 重復(fù)Step2-Step3直至收斂, 并計(jì)算λn對應(yīng)的BIC值, 其表達(dá)式為
(11)
其中,df表示λn對應(yīng)模型所選擇的變量個數(shù).
注定理1表明所提出的方法可以選出正確的模型, 同時實(shí)現(xiàn)對重要變量回歸系數(shù)的參數(shù)估計(jì), 稱為Oracle性質(zhì)[2].
定理1的證明:
(12)
成立即可. 根據(jù)表達(dá)式, 有
(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分為兩部分計(jì)算, 有
(20)
記
(21)
其中由(A3)知
又因?yàn)?/p>
(24)
(25)
其中
又
(28)
I2=-Op(pn)‖u‖2
(29)
由(19),(29)式可得, (14)式的值由(29)式控制, 小于0. 易知(13)式中的第二項(xiàng)以nαn2‖u‖ +nbnαn2‖u‖2為界, 因此可以找到一個足夠大的D, 使得(13)式的值完全由(29)式?jīng)Q定. (12)式得證.
(30)
(31)
由(A7)可知, (31)式的符號完全由βj的符號決定. (30)式得證.
即
(32)
(33)
定理證畢.
為了研究所提方法的有限樣本性質(zhì), 本文比較了不同的懲罰方法及相關(guān)結(jié)構(gòu)下所提出方法的效果. 數(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的等相關(guān)結(jié)構(gòu)矩陣.
情形2k=2,mi服從參數(shù)為(3, 6)的均勻分布,Ri是參數(shù)為0.9的AR(1)結(jié)構(gòu)矩陣. 其余設(shè)置和情形1一樣.
情形3pn=30,n=100, 200.βn=(-3, 5, 0, 0, 4, 0, 0, 2, 0, 0, …, 0)T. 其余設(shè)置和情形1一樣.
表1 情形1模擬結(jié)果
表2 情形2模擬結(jié)果
表3 情形3模擬結(jié)果
(i) SCAD和MCP兩種懲罰方法并無明顯的優(yōu)劣之分. FN均為0, 表示所有重要的變量都被識別, FP接近0, 表明噪音變量被選擇的可能性很小;
(ii) 在情形1和情形3中,τ=0.9時, Prob等于1, 而τ=0.5時, Prob的值接近0. 這表明所提出的估計(jì)量PGEEE可以在不同的τ水平下, 有效識別出正確的模型, 刻畫數(shù)據(jù)中的異方差結(jié)構(gòu);
(iii) 在不同的τ水平下, 即使選擇的變量相同, 參數(shù)估計(jì)值也可能不同(情形2). 在此情形下, 估計(jì)量的MSE和MAE隨著樣本量增大而減小, 表示該方法可以在識別出異方差的同時實(shí)現(xiàn)回歸參數(shù)的一致估計(jì);
(iv) 對比情形1和情形3, 協(xié)變量維數(shù)pn從10增加至30, 結(jié)果顯示模型中噪音變量數(shù)量增加時, PGEEE估計(jì)表現(xiàn)依然較好, 且估計(jì)量MSE減小, 表明該方法可以用于分析高維數(shù)據(jù), 排除無關(guān)變量, 識別出重要變量.
(v) 考慮相關(guān)結(jié)構(gòu)時估計(jì)量的表現(xiàn)總體上優(yōu)于獨(dú)立(IND)的情形. 即使相關(guān)結(jié)構(gòu)被誤判后, 參數(shù)估計(jì)效果依然很好, 尤其使用UN結(jié)構(gòu)時.
數(shù)據(jù)來自1976年至1982年間對美國經(jīng)濟(jì)收入動態(tài)的面板研究, 包含了連續(xù)7年595名民眾的工資水平, 屬于平衡數(shù)據(jù), 更多詳細(xì)信息參考文獻(xiàn)[15]. 該研究中, 協(xié)變量包括工作經(jīng)歷E, 工作時間W, 工作職業(yè)O(藍(lán)領(lǐng)取1, 否則0), 工作行業(yè)I(制造業(yè)取1, 否則0), 居住地S(居住在南部取1, 否則0), 種族B(黑人取1, 否則0), 是否住在都市統(tǒng)計(jì)區(qū)A(如果是取1, 否則0), 是否結(jié)婚M(結(jié)婚取1, 否則0), 性別F(女性取1, 否則0), 勞動保障U(簽合同取1, 否則0) 及受教育程度D, 響應(yīng)變量為對數(shù)變換后的工資水平.
表4給出了τ=0.01,0.5,0.95下參數(shù)的PGEEE估計(jì), 其中τ=0. 5對應(yīng)經(jīng)典的均值回歸估計(jì). 分析結(jié)果可知, 不同的懲罰方法和不同的相關(guān)結(jié)構(gòu)選出的變量基本一致. 可以看到, 在3個水平下均被選擇的變量有O,B,F,D; 均未被選擇的變量有W. 截距項(xiàng),B,F的系數(shù)估計(jì)隨著τ不同而變化, 圖1a,b為不同種族及性別對應(yīng)的工資隨時間變化的箱線圖. 男性的工資明顯高于女性, 白人的工資明顯高于黑人. 在τ=0.01時,E被認(rèn)為是噪音變量, 而在τ=0.5和0.95時被認(rèn)為是重要變量. 在τ=0.95時, 除了獨(dú)立結(jié)構(gòu)下MCP估計(jì)外, 工作行業(yè)I, 居住地S, 是否結(jié)婚M, 勞動保障U均被剔除在模型外; 而在τ=0.01 和0.5時則被認(rèn)為是重要變量. 圖1c,d,e,f為這些變量對應(yīng)的工資分布箱線圖. 以變量S為例, 可以看到, 在低分位點(diǎn)時, 居住在北部的工資要明顯高于南部, 但是在高分位點(diǎn)時, 兩者的區(qū)別并不明顯, 這與PGEEE的估計(jì)結(jié)果相吻合. 由此可見, 該方法比采用普通最小二乘估計(jì)(τ=0.5)挖掘出了更多的信息.
圖1 工資箱線圖
表4 工資數(shù)據(jù)參數(shù)估計(jì)結(jié)果
本文基于expectile提出了高維縱向數(shù)據(jù)的PGEEE估計(jì)量, 在實(shí)現(xiàn)模型變量選擇的同時, 對模型的回歸系數(shù)進(jìn)行估計(jì). 在正則條件下本文建立了PGEEE估計(jì)量的Oracle性質(zhì). 數(shù)值模擬結(jié)果顯示, MCP與SCAD懲罰及不同的協(xié)方差結(jié)構(gòu)在變量選擇方面并無明顯差異. 相較于獨(dú)立結(jié)構(gòu), 考慮相關(guān)結(jié)構(gòu)時回歸系數(shù)的估計(jì)效率更高. 多數(shù)情況下, 不確定結(jié)構(gòu)(UN)的PGEEE估計(jì)量具有較好的估計(jì)精度. 最后建立工資數(shù)據(jù)的PGEEE模型, 可以看到在不同的τ水平下, 影響工資的因素有所區(qū)別, 同一個因素影響程度也可能不同. 這表明PGEEE可以有效識別數(shù)據(jù)中的異質(zhì)結(jié)構(gòu), 比經(jīng)典的懲罰估計(jì)方程估計(jì)(PGEE)挖掘出更豐富的信息, 更合理地分析了工資的影響因素.