盛利利,亢芳圓

(北京信息科技大學 理學院,北京 100192)

0 引言

縱向數(shù)據(jù)應用廣泛,常見于醫(yī)學、經(jīng)濟學、金融學等領(lǐng)域。在許多學科,如經(jīng)濟學中,縱向數(shù)據(jù)又叫做面板數(shù)據(jù)。在縱向數(shù)據(jù)中,同一個體的不同觀測點的數(shù)據(jù)是相關(guān)的,而來自不同個體之間的觀測是獨立的。因此,在估計這類數(shù)據(jù)時,考慮具有相關(guān)性的回歸方法非常有必要。

針對縱向數(shù)據(jù)的特點,許多統(tǒng)計學者用參數(shù)模型進行了研究,如Nelder等[1]提出廣義線性模型,Liang等[2]提出廣義估計方程 (generalized estimating equation,GEE) 方法。參數(shù)回歸模型的估計理論目前已相對成熟,但是模型只考慮參數(shù)部分,在處理實際的問題時,適應性和解釋能力不足,所以在模型中加入非線性部分是一種增強模型適應性的可行方法。如Engle等[3]提出的部分線性模型,該模型結(jié)合了參數(shù)模型和非參模型的優(yōu)點,受到了不同領(lǐng)域?qū)W者的廣泛關(guān)注。

對于縱向數(shù)據(jù)下的部分線性模型,Zeger等[4]首次進行研究,并將其應用于HIV陽轉(zhuǎn)者CD4細胞數(shù)量的時間演變數(shù)據(jù)。Lind等[5]使用 profile 核估計法估計興趣參數(shù)。He等[6]對非參部分進行了B 樣條估計,提出了M估計法。Huang等[7]利用樣條逼近和GEE對其參數(shù)部分和非參數(shù)部分進行估計。薛留根等[8]利用經(jīng)驗似然方法構(gòu)造了參數(shù)分量與非參數(shù)分量的置信區(qū)間。Bai等[9]將二次推斷函數(shù) (quadratic inference function,QIF) 方法應用到縱向數(shù)據(jù)部分線性模型的統(tǒng)計推斷。文獻[4-9]的研究都涉及工作協(xié)方差矩陣的求解,但是真實的工作協(xié)方差矩陣是未知的,這就需要額外的方法來估計它。

為了簡化估計中需要的參數(shù)和工作協(xié)方差矩陣結(jié)構(gòu),Yao等[10]提出了一種在縱向數(shù)據(jù)背景下,基于Cholesky分解和profile最小二乘法技術(shù)估計非參數(shù)模型的方法,該方法避免了估計工作協(xié)方差矩陣,并且證明了所提出的估計是漸進有效的,就如真實工作協(xié)方差矩陣是先驗已知的一樣;隨后牟婷[11]將該方法推廣至縱向數(shù)據(jù)下的部分線性模型,同Yao等[10]提出的估計過程一樣運用了Cholesky分解和profile最小二乘法技術(shù),使用局部多項式估計方法處理非參數(shù)部分。田瑞琴等[12]基于Cholesky分解,提出了一種縱向數(shù)據(jù)下的半?yún)?shù)聯(lián)合均值協(xié)方差模型貝葉斯估計,運用B樣條來逼近非參數(shù)部分。文獻[10-12]巧妙利用Cholesky分解技術(shù)處理縱向數(shù)據(jù)具有的組內(nèi)相關(guān)性,避免了估計工作協(xié)方差矩陣結(jié)構(gòu),但是大都采用樣條函數(shù)處理非參部分并結(jié)合較復雜的迭代算法,才能得到興趣參數(shù)的估計值。

為了解決這個問題,可將最小二乘核估計方法應用于部分線性模型的估計中,可以不進行迭代算法,得到興趣參數(shù)的顯示解,計算簡便。如最早由Denby[13]提出的獨立數(shù)據(jù)下的部分線性模型的最小二乘核估計方法,運用了核估計來處理非參數(shù)部分,將非參數(shù)部分線性近似表示,最后運用最小二乘方法估計興趣參數(shù)。田萍等[14]將最小二乘核估計方法運用到縱向數(shù)據(jù)部分線性模型中,提出了廣義的最小二乘核估計方法。

本文結(jié)合最小二乘核估計方法和Cholesky分解技術(shù),提出一種縱向數(shù)據(jù)下部分線性模型的估計方法,用最小二乘核估計處理參數(shù)和非參數(shù)部分的估計,用Cholesky分解處理縱向數(shù)據(jù)組內(nèi)相關(guān)性,得到興趣參數(shù)的顯式解。數(shù)值模擬表明本文方法在估計參數(shù)部分和非參數(shù)部分具有良好的估計性質(zhì),并將其應用于空氣質(zhì)量數(shù)據(jù),表明估計方法在實際應用也有較好的可行性。

1 模型的提出及估計

假設(shè)有n個觀測個體,對每個個體有J次觀測,則有觀測值(Xij,Yij,Tij)(i=1,2,…,n;j=1,2,…,J)。其中Xij、Yij為第i個個體在時間為Tij時的協(xié)變量和因變量,一類重要的部分線性模型有如下形式:

Yij=βTXij+g(Tij)+εij

(1)

式中:β=(β1,β2,…,βp)T;Tij為一維協(xié)變量;g(·)為未知平滑函數(shù);εij為隨機誤差。記第i個個體的誤差向量為εi=(εi1,εi2,…,εiJ)T,{εi}相互獨立,均值和方差分別滿足E(εi)≈0,Var(εi)=C。

為了處理模型(1)中個體內(nèi)部的隨機誤差協(xié)方差矩陣C難以求解的問題,使用Yao等[10]提出的方法,用Cholesky分解將隨機誤差的協(xié)方差矩陣C轉(zhuǎn)變?yōu)閷蔷仃?進而不再估計協(xié)方差矩陣C。這里假設(shè)cov(εi|Xi,Ti)=C,存在一個對角線元素均為1的下三角矩陣Φ,滿足

cov(Φεi)=Φcov(εi)ΦT=D

(2)

εi1=ei1

εij=Φj,1εi,1+Φj,2εi,2+…+Φj,j-1εi,j-1+eij

i=1,2,…,n;j=2,3,…,J

假設(shè){ε1,ε2,…,εn}已知,得到帶有不相關(guān)誤差eij的部分線性模型:

Yi1=g(Ti1)+βTXi1+ei1

Yij=g(Tij)+βTXij+Φj,1εi,1+Φj,2εi,2+…+

Φj,j-1εi,j-1+eij

i=1,2,…,n;j=2,3,…,J

(3)

Yi1=g(Ti1)+βTXi1+ei1

i=1,2,…,n;j=2,3,…,J

令

(4)

式中:

(5)

至此,通過Cholesky分解技術(shù)已將具有組內(nèi)相關(guān)性的模型(1)轉(zhuǎn)變?yōu)榫哂薪M內(nèi)獨立性的模型(5)。

為了得到模型(5)中未知的G(T)和φ,利用最小二乘核估計方法估計模型(5)的參數(shù)和非參數(shù)部分。

假設(shè)φ已知,令

式中:K(·)為核函數(shù);h為窗寬。

其中

有

(6)

則估計的誤差方程為

H-1(I-W)Y

(7)

其中

(8)

(9)

2 數(shù)值模擬

表1 不同回歸函數(shù)下的參數(shù)估計Table 1 estimators under different regression functions

為了體現(xiàn)平滑函數(shù)g(t)擬合的優(yōu)良特性,在平滑函數(shù)g(t)分別為cos(2πt)、e-4t和個體數(shù)為30時,將g(t)估計曲線和真實曲線作比較,見圖1。

圖1 g(t)的真實曲線和估計曲線Fig.1 The real curve and estimated curve of g(t)

從圖1可以看出,真實曲線與估計曲線基本重合,擬合有很好的效果。

綜上,本文的估計方法在有限樣本模擬中有較好的表現(xiàn),說明本文的估計方法是可行的。

3 空氣質(zhì)量數(shù)據(jù)的應用

本節(jié)將本文提出的方法應用到來自 “空氣質(zhì)量在線監(jiān)測分析平臺”的空氣質(zhì)量數(shù)據(jù),數(shù)據(jù)選取北京2014年1月到2021年12月的空氣質(zhì)量數(shù)據(jù),數(shù)據(jù)包含:PM2.5質(zhì)量濃度、PM10質(zhì)量濃度、NO2質(zhì)量濃度、CO質(zhì)量濃度、SO2質(zhì)量濃度、O3質(zhì)量濃度的月平均數(shù)據(jù)。數(shù)據(jù)來源鏈接為:https://www.aqistudy.cn/historydata/。

通常一個地區(qū)的PM2.5的形成與當?shù)匾欢螘r間內(nèi)的PM10質(zhì)量濃度、NO2質(zhì)量濃度、CO質(zhì)量濃度、SO2質(zhì)量濃度、O3質(zhì)量濃度有關(guān)。在應用中,將PM2.5質(zhì)量濃度作為模型(1)中的響應變量Y,將其余氣體的質(zhì)量濃度作為模型(1)中的協(xié)變量X或T。為了確定參數(shù)部分的協(xié)變量X和非參數(shù)的協(xié)變量T,分別作PM2.5質(zhì)量濃度與PM10質(zhì)量濃度、NO2質(zhì)量濃度、CO質(zhì)量濃度、SO2質(zhì)量濃度、O3質(zhì)量濃度的散點圖,如圖2～6所示。

圖2 PM2.5與PM10質(zhì)量濃度散點圖Fig.2 Scatter plot of mass concentration of PM2.5 and PM10

圖3 PM2.5與NO2質(zhì)量濃度散點圖Fig.3 Scatter plot of mass concentration of PM2.5 and NO2

圖4 PM2.5與CO質(zhì)量濃度散點圖Fig.4 Scatter plot of mass concentration of PM2.5 and CO

圖5 PM2.5與SO2質(zhì)量濃度散點圖Fig.5 Scatter plot of mass concentration of PM2.5 and SO2

圖6 PM2.5與O3質(zhì)量濃度散點圖Fig.6 Scatter plot of mass concentration of PM2.5 and O3

觀察圖2～6發(fā)現(xiàn),PM2.5質(zhì)量濃度分別與PM10質(zhì)量濃度、NO2質(zhì)量濃度、CO質(zhì)量濃度、SO2質(zhì)量濃度有明顯的線性關(guān)系,而與O3質(zhì)量濃度沒有明顯的線性關(guān)系。因此,將PM10質(zhì)量濃度、NO2質(zhì)量濃度、CO質(zhì)量濃度、SO2質(zhì)量濃度作為參數(shù)部分的協(xié)變量,O3質(zhì)量濃度作為非參數(shù)部分的協(xié)變量。

將北京每一年的每一個月的空氣質(zhì)量數(shù)據(jù)作為一個觀測點,假定每年中月與月之間具有相關(guān)性,年與年之間相互獨立,適用于本文提出的縱向數(shù)據(jù)估計方法。

為了體現(xiàn)本文方法與其它常用方法的實際應用差異,將本文方法、半?yún)?shù)回歸(不考慮組內(nèi)相關(guān)性的部分線性模型方法)、參數(shù)回歸(不考慮組內(nèi)相關(guān)性的線性模型方法),分別應用到北京市空氣質(zhì)量數(shù)據(jù)中,得到三種模型關(guān)于真實PM2.5估計的折線如圖7所示,同時計算每種方法的PM2.5估計值與真實值之間的均方誤差,見表2。

圖7 真實數(shù)據(jù)與三種估計的折線圖Fig.7 Real data and three estimated line charts

表2 不同回歸方法關(guān)于PM2.5估計值與真實值的均方誤差Table 2 The mean square error of PM2.5 estimated value and real value by different regression methods

從圖7可看出,本文方法和半?yún)?shù)回歸方法擬合的折線與真實PM2.5折線無明顯差異,而參數(shù)方法擬合的折線與真實PM2.5折線差異較大,說明參數(shù)模型不適用擬合本節(jié)的空氣質(zhì)量數(shù)據(jù)。

從表2可看出,本文方法均方誤差最小,說明如果忽略數(shù)據(jù)每一年中月與月之間的相關(guān)性,會對擬合有較大的影響,簡單地考慮數(shù)據(jù)之間的關(guān)系對擬合的效果有很大的折損。本文方法回歸擬合中參數(shù)部分的協(xié)變量系數(shù)分別為0.143,0.538,28.571,0.206,說明PM10質(zhì)量濃度、NO2質(zhì)量濃度、CO質(zhì)量濃度、SO2質(zhì)量濃度這四種因素都影響著PM2.5的形成。

綜上可得,PM2.5的產(chǎn)生與PM10質(zhì)量濃度、NO2質(zhì)量濃度、CO質(zhì)量濃度、SO2質(zhì)量濃度有密切的關(guān)系;北京市空氣質(zhì)量數(shù)據(jù)中,PM2.5的數(shù)值逐年下降,北京的空氣環(huán)境治理有明顯的成效。將空氣質(zhì)量數(shù)據(jù)作為縱向數(shù)據(jù)并使用本文提出的估計方法,能夠得到較好的擬合效果。

4 結(jié)束語

本文將Cholesky 分解和最小二乘核估計方法相結(jié)合,給出了新估計方法的具體估計過程,用Cholesky 分解來處理縱向數(shù)據(jù)的組內(nèi)相關(guān)性,用最小二乘核估計來處理部分線性模型中的非參數(shù)部分,由于可以得到顯式解,所以該估計方法可以較為迅速地得出估計結(jié)果,比一般的迭代方法有明顯的優(yōu)勢,并且不再對縱向數(shù)據(jù)中的工作協(xié)方差矩陣進行估計,簡化的估計程序又節(jié)約了計算成本。通過數(shù)據(jù)模擬說明該估計方法的可靠性。本文將估計方法運用于空氣質(zhì)量數(shù)據(jù)的分析中,結(jié)果表明該估計方法適用于常見數(shù)據(jù)的應用分析。這也為分析其它類型的縱向數(shù)據(jù)提供了解決思路。