函數(shù)型部分線性模型的參數(shù)估計(jì)及其變量選擇

袁曉惠，曹儒雅，金宛霖

(長春工業(yè)大學(xué) 數(shù)學(xué)與統(tǒng)計(jì)學(xué)院，吉林 長春 130012)

0 引言

函數(shù)型數(shù)據(jù)最初是由J.Ramsay[1]提出的，是一種無限維的數(shù)據(jù).函數(shù)型數(shù)據(jù)是指在緊區(qū)域上的實(shí)值函數(shù)的一個(gè)實(shí)現(xiàn)[2].由離散數(shù)據(jù)擬合得到光滑曲線是函數(shù)型數(shù)據(jù)的研究對象[3].由于函數(shù)型數(shù)據(jù)的無限維屬性，使得以往的多元統(tǒng)計(jì)方法不再適用，這為函數(shù)型數(shù)據(jù)的發(fā)展提供了先決條件.此外，由于函數(shù)型數(shù)據(jù)是無窮維的，在現(xiàn)實(shí)中有一些數(shù)據(jù)是不能完全被觀測到的，只能觀測到其有限的離散時(shí)間點(diǎn).因此可將函數(shù)型數(shù)據(jù)根據(jù)數(shù)量大小分為稠密數(shù)據(jù)和稀疏數(shù)據(jù)，根據(jù)取值形式分為均衡數(shù)據(jù)和非均衡數(shù)據(jù).根據(jù)樣本觀測量的大小及分析函數(shù)型數(shù)據(jù)的工具的不同，可將研究函數(shù)型數(shù)據(jù)的學(xué)者分為三個(gè)學(xué)派，即經(jīng)典學(xué)派、隨機(jī)學(xué)派和法國學(xué)派.

近年來，關(guān)于函數(shù)型數(shù)據(jù)分析已經(jīng)積累了大量的成果.常用的函數(shù)型數(shù)據(jù)分析的工具是函數(shù)型線性模型，主要有兩種：一是標(biāo)量型響應(yīng)變量的函數(shù)型線性模型，二是函數(shù)型響應(yīng)變量的函數(shù)型線性模型.本文只考慮響應(yīng)變量為標(biāo)量型的函數(shù)型回歸模型，再根據(jù)協(xié)變量的數(shù)據(jù)類型，可以將函數(shù)型回歸模型分為兩類：協(xié)變量為函數(shù)型數(shù)據(jù)的回歸模型和協(xié)變量為函數(shù)型與標(biāo)量型混合數(shù)據(jù)的回歸模型.

比較經(jīng)典的協(xié)變量為函數(shù)型數(shù)據(jù)的函數(shù)型回歸模型是函數(shù)型線性模型.該模型是研究最早和最廣泛的模型之一.H.Cardot[4]首次研究了該模型的計(jì)算問題，并于2003年基于樣條基底展開和主成分降維的估計(jì)，推導(dǎo)了基于樣條展開估計(jì)的漸近收斂性質(zhì)；T.T.Cai等[5]在理論上推導(dǎo)了函數(shù)線性模型預(yù)測的收斂速度，在一定條件下，得出了預(yù)測能達(dá)到參數(shù)收斂速率的結(jié)論；H.G.Muller等[6]研究了基于主成分得分的函數(shù)加性模型的估計(jì)和理論性質(zhì)；F.Yao等[7]研究了函數(shù)型數(shù)據(jù)二次回歸模型；在再生核希爾伯持空間上，M.Yuan等[8]提出懲罰估計(jì)相對于主成分回歸能在較弱的條件下達(dá)到最優(yōu)收斂速率；T.T.Cai等[9]在再生核希爾伯特空間里研究了預(yù)測的最優(yōu)收效速率和自適應(yīng)性；Z.Lin等[10]利用最小二乘基于樣條技術(shù)進(jìn)行了函數(shù)型回歸模型的局部稀疏估計(jì)；Q.Wang等[11]提出了基于函數(shù)充分性降維基底來估計(jì)函數(shù)的線性模型，并推導(dǎo)了相關(guān)理論性質(zhì).

比較經(jīng)典的協(xié)變量為混合數(shù)據(jù)的函數(shù)型回歸模型是函數(shù)型部分線性模型.此模型由D.Zhang等[12]提出；H.Shin[13]主要研究的是模型基于函數(shù)型主成分的估計(jì)問題，并得到了參數(shù)部分的漸進(jìn)正態(tài)性和函數(shù)系數(shù)的收斂速度；H.Shin等[14]使用基于FPCA 的最小二乘方法和嶺回歸方法研究了該模型的估計(jì)和預(yù)測，并證明了在一定的條件下斜率函數(shù)的估計(jì)能達(dá)到最優(yōu)收斂速度；丁輝[15]對函數(shù)型部分線性單指標(biāo)模型、函數(shù)型部分線性回歸模型等幾類函數(shù)型回歸模型的估計(jì)與預(yù)測展開了研究.

在實(shí)際應(yīng)用中，遇到的協(xié)變量大部分是多維的，而在這些協(xié)變量中有一些是無關(guān)變量，變量選擇是剔除無關(guān)變量的一個(gè)很好的工具.常用的變量選擇方法有嶺回歸、Lasso、SCAD等.Lasso變量選擇方法的優(yōu)點(diǎn)是，它既改進(jìn)了傳統(tǒng)變量選擇方法在模型選擇方面的不足之處，同時(shí)又保留了一些優(yōu)良的性質(zhì)[16].此外，嶺回歸和 Lasso 回歸方法也可以消除自變量間存在的共線性問題[17]，同時(shí)具有一定的穩(wěn)定性[18].SCAD方法則改進(jìn)了Lasso方法在模型應(yīng)用上的一些局限.在函數(shù)型部分線性模型中，也可采用適當(dāng)?shù)淖兞窟x擇方法，剔除無關(guān)變量.

本文的估計(jì)方法和結(jié)論不僅豐富了函數(shù)型回歸模型的研究，給函數(shù)型回歸模型的發(fā)展提供了更多的參考，更將有助于以后在經(jīng)濟(jì)學(xué)、醫(yī)學(xué)、生物學(xué)等領(lǐng)域的發(fā)展與探索.

1 函數(shù)型部分線性模型的估計(jì)方法

1.1 參數(shù)估計(jì)

函數(shù)型部分線性模型是研究函數(shù)型數(shù)據(jù)的一個(gè)應(yīng)用較為廣泛的模型[12]，形式為

(1)

其中：Yi為響應(yīng)變量；X1(t)，X2(t)，…，Xn(t)是定義在區(qū)間[0，T]上的函數(shù)型協(xié)變量；Zi=(Zi1，…，Zip)T為p維協(xié)向量；β(t)為斜率函數(shù)；α=(α1，…，αp)T為Zi的系數(shù)；εi為誤差項(xiàng)，且εi～N(0，σ2).此模型即描述了標(biāo)量型響應(yīng)變量和函數(shù)型與標(biāo)量型的混合協(xié)變量之間的關(guān)系.

用B樣條方法表示β(t)，即

其中：B(t)=(B1(t)，B2(t)，…，BM+d(t))T是B樣條基函數(shù)；b=(b1，b2，…，bM+d)T是相應(yīng)的系數(shù).但是當(dāng)M較大時(shí)，直接使用B樣條方法近似，會使得估計(jì)值產(chǎn)生較大的波動性，因此，在目標(biāo)函數(shù)中加入粗糙度懲罰函數(shù)，即

(2)

其中Dmβ為β(t)的m階導(dǎo)數(shù)，且m

此外，當(dāng)M較大時(shí)，可以通過加入光滑消邊絕對偏離(SCAD)懲罰函數(shù)pλ1(u)來實(shí)現(xiàn)斜率函數(shù)空子區(qū)間的尋找，同時(shí)加入光滑消邊絕對偏離(SCAD)懲罰函數(shù)對α中零分量的尋找.即最小化下面的目標(biāo)函數(shù)：

其中

通過最小化上面的目標(biāo)函數(shù)，可以實(shí)現(xiàn)斜率函數(shù)空子區(qū)間與非空子區(qū)間的尋找和標(biāo)量型系數(shù)零分量的識別.為了計(jì)算方便，引入一些符號，令

令

Y=(Y1，…，Yn)T，Z=(Z1，…，Zn)T，X(t)=(X1(t)，…，Xn(t))T，

U=(uij)n×(M+d)，V=(vij)(M+d)×(M+d)，G=(U，Z)，g=(b，α).

由此可得

計(jì)算得到

目標(biāo)函數(shù)的第二項(xiàng)可轉(zhuǎn)化為γ‖Dmβ‖2=γbTVb.為了方便計(jì)算目標(biāo)函數(shù)的第三項(xiàng)，利用文獻(xiàn)[10]中的定理1，fSCAD懲罰項(xiàng)可被近似為

利用LQA方法來計(jì)算，即

其中：

但由于H(β(0))不依賴于β，因此對目標(biāo)函數(shù)沒有影響，即可寫為

目標(biāo)函數(shù)第四項(xiàng)中的pλ2(|αj|)可以進(jìn)行泰勒展開，即

(3)

對式(3)求一階導(dǎo)，有

求二階導(dǎo)，有

令

=(GTG+nSS)-1GTY.

1.2 算法步驟

為了更清楚地表達(dá)上述估計(jì)思想，列出如下估計(jì)步驟：

1.3 調(diào)節(jié)參數(shù)選擇

為了實(shí)現(xiàn)上述估計(jì)，需要進(jìn)行調(diào)節(jié)參數(shù)的選擇.經(jīng)過多次研究分析，γ=1e-7時(shí)效果最好.下面進(jìn)行λ1，λ2的選擇.

GCV準(zhǔn)則定義為

其中e(λ1)=tr[PG{(λ1)}]，PG{(λ1)}=GGT.

BIC準(zhǔn)則定義為

2 模擬分析

(Xi(t)，Zi，Yi)由如下模型產(chǎn)生：

第一種β(t)=0，即只存在空子區(qū)間.

β(t)估計(jì)值的表現(xiàn)通過積分平方誤(ISE)來評價(jià)，即

其中l(wèi)0和l1表示空子區(qū)間和非空子區(qū)間的長度，N(β)和S(β)則表示對應(yīng)的空子區(qū)間和非空子區(qū)間.α估計(jì)值的表現(xiàn)則通過偏差和均方誤差來體現(xiàn).模擬均重復(fù)1 000次.

第一種斜率函數(shù)的結(jié)果如圖1所示.

圖1 情形1中β(t)的真實(shí)曲線圖及其估計(jì)曲線

由此可見，當(dāng)斜率函數(shù)為零時(shí)，Spline方法可使得估計(jì)值在零附近上下波動，而SLOS方法可使斜率函數(shù)的估計(jì)完全變?yōu)榱?因此，在只存在空子區(qū)間的情況下，SLOS方法比Spline方法效果更好.具體分析結(jié)果如表1所示：

表1 情形1估計(jì)的偏差和(均方誤差)或ISE0、ISE1和(標(biāo)準(zhǔn)差)

由表1可分析得，SLOS方法的ISE0比Spline方法的ISE0要小很多，這表明SLOS可以使斜率函數(shù)完全被懲罰為零.而對α的估計(jì)方面，SLOS方法能夠更好的識別零向量，Spline方法則不能達(dá)到這樣的效果.因此，總體而言，SLOS方法要優(yōu)于Spline方法.

第二種斜率函數(shù)的結(jié)果如圖2所示.

圖2 情形2中β(t)的真實(shí)曲線圖及其估計(jì)曲線

由此可見，當(dāng)斜率函數(shù)同時(shí)存在空子區(qū)間與非空子區(qū)間的時(shí)候，SLOS方法可以很好的識別空子區(qū)間，使其被估計(jì)為零.與Z.Lin[10]結(jié)論一致.具體分析結(jié)果如表2所示.

表2 情形2估計(jì)的偏差和(均方誤差)或ISE0、ISE1和(標(biāo)準(zhǔn)差)

由表2可分析得，SLOS方法的ISE0比Spline方法的ISE0要小很多，即SLOS方法可以很好的識別空子區(qū)間.而隨著樣本量的增加，兩種方法的ISE1之間的差值逐漸減小.而對α的估計(jì)方面，SLOS方法能夠更好地識別零向量，因此整體上，SLOS方法要比Spline方法更加適用于此模型.

第三種斜率函數(shù)的結(jié)果如圖3所示.

圖3 情形3中β(t)的真實(shí)曲線圖及其估計(jì)曲線

由此可見，當(dāng)斜率函數(shù)不為零，即只存在非空子區(qū)間時(shí)，Spline方法的估計(jì)曲線與真值曲線趨于一致，明顯優(yōu)于SLOS方法.因此，在這種情況下，SCAD懲罰函數(shù)將不再適用.具體分析結(jié)果如表3所示.

由表3可分析得，Spline方法的ISE1比SLOS方法的ISE1要小很多，即在只存在非空子區(qū)間的情況下，Spline方法更加適用.在α的估計(jì)方面，SLOS方法的效果更好.這表明，SLOS方法更加適用于存在空子區(qū)間或者零向量的情況，而只有非空子區(qū)間或者非零向量的情況下，Spline方法則更優(yōu)越.

表3 情形3估計(jì)的偏差和(均方誤差)或ISE0，ISE1和(標(biāo)準(zhǔn)差)

3 實(shí)例分析

本節(jié)將利用本文中所述的模型和方法對2018年中國主要的31個(gè)城市的空氣質(zhì)量進(jìn)行研究分析，該數(shù)據(jù)來源于https：//quotsoft.net/air/以及中國統(tǒng)計(jì)年鑒.空氣質(zhì)量數(shù)據(jù)是由31個(gè)城市樣本組成，其中每個(gè)樣本包括PM2.5(μg/m3)含量，CO(mg/m3)含量，氣溫(℃)，濕度(%)，降水量(mm)，日照時(shí)數(shù)(h)，工業(yè)廢水排放量(t)，工業(yè)SO2(t)排放量，工業(yè)氨氮排放量(t).其中CO含量是由一年中每天24 h滑動均值構(gòu)成的.首先對降水量、日照時(shí)數(shù)、工業(yè)廢水排放量、工業(yè)SO2排放量和工業(yè)氨氮排放量等數(shù)據(jù)進(jìn)行標(biāo)準(zhǔn)化，擬合的模型為

其中：Y是PM2.5含量；X(t)是CO含量；Z是氣溫、濕度、降水量、日照時(shí)數(shù)、工業(yè)廢水排放量、工業(yè)SO2排放量和工業(yè)氨氮排放量.CO含量見圖4.

由圖4分析可得，在4月份到9月份之間空氣中CO含量相對穩(wěn)定，而在1—2月和11—12月之間CO含量則幾乎兩倍于4—9月份.這可以給環(huán)保部門在進(jìn)行措施制定時(shí)提供一些依據(jù).此外，β(t)估計(jì)值如圖5所示.

圖4 中國一年31個(gè)城市日均CO質(zhì)量濃度

由圖5分析可知，β(t)值整體波動趨于零，而在前端和后端有較大起伏，即在一年中秋冬季節(jié)CO含量對PM2.5含量有較大影響，此外，空氣中CO含量對PM2.5含量的影響處于一個(gè)相對穩(wěn)定的區(qū)間內(nèi).

圖5 β(t)的估計(jì)曲線

得到α估計(jì)值如表4所示.

表4 α的估計(jì)值

由表4分析可知，降水量對于PM2.5含量呈負(fù)相關(guān)，即降水量越多，PM2.5含量越少，而日照時(shí)數(shù)、工業(yè)廢水排放量和工業(yè)SO2排放量對于PM2.5含量的影響則呈正相關(guān)性.氣溫、濕度與工業(yè)氨氮排放量則相對沒有影響.

以上信息不僅對相關(guān)部門、工廠、居民等在日常生活中有針對性、目的性、計(jì)劃性和科學(xué)性地減少污氣、廢水排放，節(jié)能環(huán)保，廢物利用，降低PM2.5含量等起到?jīng)Q定性作用，也更有助于今后的研究.

4 結(jié)語

本文基于最小二乘和光滑樣條法，加入光滑消邊絕對偏離懲罰函數(shù)(SCAD)實(shí)現(xiàn)空子區(qū)間的尋找，再對標(biāo)量型系數(shù)進(jìn)行變量選擇，得到斜率函數(shù)的局部稀疏估計(jì)和標(biāo)量型系數(shù)中零分量的識別.在統(tǒng)計(jì)模擬部分，所提出的估計(jì)方法展現(xiàn)出良好的有限樣本表現(xiàn).最后的實(shí)證部分，通過對2018年中國主要的31個(gè)城市的空氣質(zhì)量進(jìn)行研究分析，再次論證了本文所述估計(jì)方法的優(yōu)良性與適用性.

此外，局部稀疏估計(jì)方法適用于很多的函數(shù)型回歸模型，而本文只研究了函數(shù)型部分線性模型.所以，局部稀疏估計(jì)方法在其他的函數(shù)型回歸模型中的應(yīng)用是未來的一個(gè)研究方向.同時(shí)，只研究了一個(gè)函數(shù)型協(xié)變量的情況和有限維的標(biāo)量型協(xié)變量的情況，而大多數(shù)情況下，函數(shù)型協(xié)變量和標(biāo)量型協(xié)變量的維數(shù)很高，甚至是無限維的.因此，高維的函數(shù)型部分線性模型的局部稀疏估計(jì)與變量選擇相結(jié)合是未來的又一個(gè)研究方向.