基于指數(shù)平方損失的縱向多折點(diǎn)回歸模型的穩(wěn)健估計(jì)與統(tǒng)計(jì)推斷①

唐銘， 李婷婷

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

折點(diǎn)回歸模型主要研究響應(yīng)變量與協(xié)變量間的分階段連續(xù)變化特征, 在金融、 經(jīng)濟(jì)、 工業(yè)、 醫(yī)學(xué)等領(lǐng)域有著重要應(yīng)用． 文獻(xiàn)[1]首次提出單折點(diǎn)回歸模型, 并基于似然函數(shù)提出折點(diǎn)參數(shù)估計(jì)的網(wǎng)格搜索法． 文獻(xiàn)[2]基于累積和統(tǒng)計(jì)量及似然比統(tǒng)計(jì)量提出折點(diǎn)效應(yīng)的假設(shè)檢驗(yàn)方法． 值得一提的是, 文獻(xiàn)[1]所提的網(wǎng)格搜索法雖然可以生成合理的估計(jì), 但計(jì)算成本卻很高． 文獻(xiàn)[3]基于泰勒展開(kāi)提出了一種局部線性平滑的參數(shù)估計(jì)方法, 在保證估計(jì)準(zhǔn)確性的同時(shí)極大提高計(jì)算效率． 文獻(xiàn)[4]首次將單折點(diǎn)回歸模型拓展到多折點(diǎn)回歸模型, 基于局部線性平滑法提出了分位數(shù)損失函數(shù)下的參數(shù)估計(jì)、 變點(diǎn)存在性檢驗(yàn)統(tǒng)計(jì)量, 及折點(diǎn)個(gè)數(shù)確定的貝葉斯信息準(zhǔn)則, 并研究了所提估計(jì)及統(tǒng)計(jì)量的大樣本性質(zhì)．

以上文獻(xiàn)的研究大多討論的是獨(dú)立同分布數(shù)據(jù)的折點(diǎn)模型． 然而, 隨著應(yīng)用領(lǐng)域的不斷拓展, 所處理的數(shù)據(jù)類型越來(lái)越復(fù)雜, 縱向數(shù)據(jù)便是復(fù)雜的數(shù)據(jù)類型之一． 針對(duì)縱向折點(diǎn)回歸模型, 文獻(xiàn)[5]在獨(dú)立工作矩陣下考慮了縱向單折點(diǎn)分位數(shù)回歸模型的估計(jì)與檢驗(yàn); 文獻(xiàn)[6]考慮了縱向數(shù)據(jù)的多折點(diǎn)分位數(shù)回歸模型． 為融合縱向數(shù)據(jù)同一個(gè)個(gè)體內(nèi)部的相關(guān)性, 文獻(xiàn)[6]基于文獻(xiàn)[7]提出的二次推斷函數(shù)方法(quadratic inference function, QIF)研究了相關(guān)結(jié)構(gòu)下縱向多折點(diǎn)分位數(shù)回歸模型的估計(jì)與統(tǒng)計(jì)推斷, 但其所能刻畫的仍是等相關(guān)和AR(1)等特殊結(jié)構(gòu)的矩陣． 文獻(xiàn)[8]提出修正的Cholesky分解方法, 該方法不局限于特殊相關(guān)結(jié)構(gòu), 且能保證估計(jì)的協(xié)方差陣的正定性, 具有更廣泛的適用性．

眾所周知, 基于經(jīng)典平方損失的估計(jì)對(duì)異常值非常敏感． 為處理包含大量異常值的數(shù)據(jù), 眾多穩(wěn)健估計(jì)的方法被相繼提出, 如Huber損失函數(shù)法[9]、 秩回歸[10]及分位數(shù)回歸方法[11]等． 文獻(xiàn)[12]提出一種新的穩(wěn)健估計(jì)方法, 即基于指數(shù)平方損失函數(shù)的參數(shù)估計(jì)方法． 該方法的顯著特征是引入一個(gè)額外的調(diào)諧參數(shù), 可通過(guò)選擇合適的調(diào)諧參數(shù)實(shí)現(xiàn)模型參數(shù)的自適應(yīng)穩(wěn)健估計(jì)． 文獻(xiàn)[13-15]關(guān)于指數(shù)平方損失函數(shù)的相關(guān)研究均表明, 基于該損失的參數(shù)估計(jì)相對(duì)于經(jīng)典的穩(wěn)健方法有著更好的表現(xiàn), 能夠獲得更好的魯棒性和有效性．

本文基于指數(shù)平方損失和修正的Cholesky分解方法研究縱向多折點(diǎn)回歸模型的參數(shù)估計(jì)及統(tǒng)計(jì)推斷．

1 模型與估計(jì)算法

1.1 縱向多折點(diǎn)回歸模型

考慮折點(diǎn)個(gè)數(shù)K及折點(diǎn)位置τ=(τ1, …,τK)T均未知的縱向多折點(diǎn)回歸模型:

(1)

其中:n為個(gè)體數(shù),mi為第i個(gè)個(gè)體的重復(fù)觀測(cè)次數(shù), (a)+=a·I(a≥0),I為示性函數(shù),Yij為響應(yīng)變量觀測(cè)值,Zij為p維協(xié)變量觀測(cè)值,Xij為有界門限變量,eij為隨機(jī)誤差, 記ei=(ei1, …,eimi)T． 由模型(1)可知, 變量Xij的回歸系數(shù)在折點(diǎn)τ=(τ1, …,τK)T處會(huì)發(fā)生變化． 記b=(b1, …,bK)T, ?= (a0,a1,bT,βT)T, 以及參數(shù)向量θ= (?T,τT)T, 其中?∈G?R2+K+p,τ∈T?DK．G,T均為緊集,D為Xij的支撐集． 下面, 首先對(duì)假設(shè)折點(diǎn)個(gè)數(shù)K已知后的參數(shù)估計(jì)算法進(jìn)行介紹, 隨后給出用于確定折點(diǎn)個(gè)數(shù)K的模型選擇方法．

1.2 基于指數(shù)平方損失的縱向多折點(diǎn)回歸模型的參數(shù)估計(jì)算法

(2)

進(jìn)而有如下的近似回歸模型

(3)

(4)

(5)

(6)

其中Φi是主對(duì)角線元素全為1的下三角矩陣, 第(j,k)個(gè)元素是如下自回歸方程系數(shù)φijk,γ的相反數(shù)

(7)

(8)

其中

本文使用Newton-Raphson迭代算法求解方程組(8)以保證估計(jì)的精度． 下面給出在指定調(diào)諧參數(shù)γ下, 基于指數(shù)平方損失的縱向多折點(diǎn)回歸模型的參數(shù)估計(jì)算法:

步驟1給定初始折點(diǎn)位置τ(0), 計(jì)算模型(3)的普通最小二乘估計(jì)作為η的初始值η(0)．

步驟2基于τ(s)及η(s), 應(yīng)用修正的Cholesky分解法估計(jì)協(xié)方差陣Vi, 使用Newton-Raphson算法迭代求解方程組(8)獲得η(s+1), 具體步驟如下:

步驟2.2.3根據(jù)(6)式及(7)式, 可得Vi, r+1． 通過(guò)下式迭代至收斂得到η(s+1, r+1):

步驟2.3重復(fù)步驟2.2直至參數(shù)收斂, 可得η(s+1)．

步驟3更新折點(diǎn)位置τ(s+1)． 通過(guò)

(9)

步驟4重復(fù)步驟2-步驟3直至參數(shù)收斂．

注1當(dāng)Vi為單位陣時(shí), 方程(5)和最小化目標(biāo)函數(shù)(4)的解等價(jià)． 因此, 在上述步驟中省略步驟2.2.1及步驟2.2.2, 設(shè)置步驟2.2.3中Vi, r+1為單位陣, 即可獲得獨(dú)立工作矩陣情況的參數(shù)估計(jì)．

1.3 最佳調(diào)諧參數(shù)γopt的選擇

1.4 折點(diǎn)個(gè)數(shù)K的確定

1.2節(jié)給出了當(dāng)折點(diǎn)個(gè)數(shù)給定的時(shí)候模型參數(shù)的估計(jì)算法． 然而實(shí)際問(wèn)題中, 折點(diǎn)個(gè)數(shù)真值K0通常是未知的． 參考文獻(xiàn)[4], 提出如下基于指數(shù)平方損失的貝葉斯準(zhǔn)則以確定折點(diǎn)個(gè)數(shù):

(10)

2 大樣本性質(zhì)

(11)

給出如下條件:

(ⅰ) E(ψγ(eij))=0, E(ψ′γ(eij))>0, 且對(duì)于任意γ>0, E(ψγ(eij)2)有界．

(ⅱ) 折點(diǎn)的真值K0及協(xié)變量Zij的維數(shù)p是固定的, 個(gè)體的觀測(cè)次數(shù)mi一致有界, 且個(gè)體數(shù)n趨于無(wú)窮．

(ⅵ) 參數(shù)空間Θ∈R2+p+2K是緊集, 參數(shù)真值η0是參數(shù)空間Θ的內(nèi)點(diǎn),Sn(η)在η0處唯一取得全局最小值．

(12)

其中

(13)

(14)

3 折點(diǎn)的存在性檢驗(yàn)

對(duì)于多折點(diǎn)回歸模型的折點(diǎn)存在性檢驗(yàn), 文獻(xiàn)[4]基于無(wú)折點(diǎn)的原假設(shè)H0:bk=0,k=1, …,K提出CUSUM(cumulative summation)型統(tǒng)計(jì)量． 本文沿用文獻(xiàn)[4]所提方法, 提出基于指數(shù)平方損失的縱向折點(diǎn)回歸模型的折點(diǎn)存在性檢驗(yàn)方法, 具體步驟如下:

步驟4重復(fù)步驟3L次, 計(jì)算

4 模擬研究

根據(jù)如下模型生成數(shù)據(jù)

(15)

情形2 考慮門限變量不為時(shí)間, 且數(shù)據(jù)存在異常值的情況． 設(shè)置Xij～U(-5, 5),tij={1, 2, …, 10}, 隨后每個(gè)觀測(cè)以20%的概率缺失以生成不平衡數(shù)據(jù)．ei服從自由度為3, 協(xié)方差陣為情形1中Σi的多元t分布, 隨后隨機(jī)選取10%的eij服從標(biāo)準(zhǔn)柯西分布．

每種情形均考慮3種折線效應(yīng), 對(duì)于情形1, 設(shè)定: ①K=1,τ=7, ?=(-2, 1, -3, 1)T; ②K=2,τ=(5, 7)T, ?=(2, 1, -4, 5, 1)T; ③K=3,τ=(2, 5, 7)T, ?=(0 , 1 , -3 , 2 , -6 , 1)T． 對(duì)于情形2、 3, 系數(shù)?設(shè)置與情形1相同, 折點(diǎn)個(gè)數(shù)和位置設(shè)置為: ①K=1,τ=-1; ②K=2,τ=(-1, 1)T; ③K=3,τ=(-4 , -1 , 4)T． 設(shè)定樣本量n=400, 重復(fù)模擬次數(shù)為100次． 使用1.2節(jié)所提算法求解最小化目標(biāo)函數(shù)(4)及方程組(8), 分別記為“ESL．IND”和“ESL．CHO”, 即獨(dú)立工作矩陣結(jié)構(gòu)及基于Cholesky分解的模型參數(shù)的估計(jì)方法．

4.1 最優(yōu)調(diào)諧參數(shù)γopt的選擇

表1給出3種情形下所有模擬中選擇的最優(yōu)參數(shù)γopt的均值． 可以看到, 無(wú)異常值的情形下,γopt均值較大; 情形2,3的γopt較小, 以降低異常值產(chǎn)生的影響． 這說(shuō)明在指數(shù)平方損失函數(shù)中, 可以通過(guò)選擇合適的調(diào)諧參數(shù)實(shí)現(xiàn)回歸系數(shù)的自適應(yīng)穩(wěn)健估計(jì)．

表1 最優(yōu)調(diào)諧參數(shù)γopt的均值

4.2 選擇相合性及參數(shù)估計(jì)

表2給出了真實(shí)折點(diǎn)個(gè)數(shù)為2時(shí), 不同Cn在3種情形下的折點(diǎn)個(gè)數(shù)的正確選擇率． 可以看到, 所有情形下, 本文所提出的折點(diǎn)個(gè)數(shù)選擇方法均具有較高的正確選擇率, 并且“ESL.CHO”能夠?qū)崿F(xiàn)更高的選擇正確率． 而且樣本量不變時(shí), 較大的Cn能略微提高正確選擇折點(diǎn)數(shù)的概率．

表2 K=2, Cn=1, log(log n), log n, 2log n時(shí)的折點(diǎn)個(gè)數(shù)正確選擇率

為說(shuō)明所提方法優(yōu)勢(shì), 與以下估計(jì)量進(jìn)行比較:

1) 方程(5)的工作協(xié)方差矩陣Vi為等相關(guān)及AR(1)的特定結(jié)構(gòu), 利用文獻(xiàn)[7]提出的QIF方法對(duì)(5)式進(jìn)行求解, 同時(shí)分別考慮指數(shù)平方損失及經(jīng)典的平方損失兩種損失函數(shù), 記所得估計(jì)量分別為“ESL．EXCH”“OLS．EXCH”“ESL．AR1”“OLS．AR1”;

2) 文獻(xiàn)[4]提出的基于分位數(shù)回歸的縱向多折點(diǎn)模型的估計(jì)量, 指定分位水平為0.5, 記為“MKQR”, 使用R程序包MultiKink實(shí)現(xiàn);

3) 文獻(xiàn)[3]提出的多折點(diǎn)模型的最小二乘估計(jì)量, 記為“SEG”, 使用R程序包segmented實(shí)現(xiàn)．

表3展示了情形2在K=2時(shí)的參數(shù)估計(jì)結(jié)果, 匯總了100次模擬中估計(jì)量的平均偏差、 標(biāo)準(zhǔn)差及均方誤差． 其余情形的參數(shù)估計(jì)的結(jié)果已上傳至Github(https: //github．com/Tangming-hub)． 可以發(fā)現(xiàn):

表3 情形2, K=2的模擬結(jié)果

1) 當(dāng)殘差向量服從正態(tài)分布時(shí)(情形1), 幾種方法所得到的估計(jì)量的估計(jì)效果相近． 然而當(dāng)數(shù)據(jù)中存在異常值時(shí)(情形2和情形3), 基于指數(shù)平方損失和分位數(shù)回歸的估計(jì)量“ESL．IND”“ESL．CHO”“ESL．EXCH”“ESL．AR1”和“MKQR”均能提供回歸參數(shù)的有效估計(jì), 其中相對(duì)于分位數(shù)回歸方法, 基于指數(shù)平方損失函數(shù)的估計(jì)量表現(xiàn)更佳, 而基于經(jīng)典的平方損失函數(shù)的估計(jì)量的估計(jì)效果均較差;

2) 僅考慮指數(shù)平方損失函數(shù)的估計(jì)方法時(shí), 相對(duì)于獨(dú)立工作矩陣的估計(jì)量“ESL．IND”, 融合組內(nèi)相關(guān)性的估計(jì)量“ESL．CHO”“ESL．EXCH”“ESL．AR1”均具有更優(yōu)良的估計(jì)效果, 且本文所提出的“ESL．CHO”在各指標(biāo)上表現(xiàn)最佳． 這說(shuō)明有效考慮縱向數(shù)據(jù)個(gè)體重復(fù)觀測(cè)有利于提高回歸參數(shù)的估計(jì)效率． 綜上, 本文所提出的估計(jì)方法, 可以為縱向多折點(diǎn)回歸模型的參數(shù)提供更為穩(wěn)健的、 有效的估計(jì)．

此外, 不同情形下K=2時(shí), 3種估計(jì)方法“ESL．CHO”“ESL．IND”“MKQR”的標(biāo)準(zhǔn)差、 標(biāo)準(zhǔn)誤、 95%的 Wald型置信區(qū)間的平均長(zhǎng)度及經(jīng)驗(yàn)覆蓋率亦上傳至Github． 由結(jié)果可見(jiàn), 3種估計(jì)量的標(biāo)準(zhǔn)差與經(jīng)驗(yàn)標(biāo)準(zhǔn)誤差接近, 并且所構(gòu)造的置信區(qū)間的經(jīng)驗(yàn)覆蓋率均在置信水平95%左右． 并且, 相對(duì)于“MKQR”方法, 本文所提方法的置信區(qū)間的長(zhǎng)度更短．

4.3 折線效應(yīng)存在性檢驗(yàn)的功效性分析

為研究本文第3節(jié)所提出的折線效應(yīng)存在性檢驗(yàn)統(tǒng)計(jì)量Tn的有限樣本性質(zhì), 考慮折點(diǎn)個(gè)數(shù)為2時(shí)檢驗(yàn)統(tǒng)計(jì)量的功效． 對(duì)于情形1, 設(shè)置-b1=b2=0,0.1,0.2,0.3; 對(duì)于情形2和情形3, 設(shè)置-b1=b2=0,0.15,0.3,0.45, 其中當(dāng)系數(shù)b1和b2都等于0時(shí), 不存在折線效應(yīng)． 設(shè)置L=300, 顯著性水平α=0.05． 表4展示所得檢驗(yàn)統(tǒng)計(jì)量均值(Mean ofTn, 簡(jiǎn)記為Mean-Tn)及經(jīng)驗(yàn)p值(Power)． 根據(jù)統(tǒng)計(jì)量的定義, 當(dāng)原假設(shè)成立時(shí), 統(tǒng)計(jì)量Tn應(yīng)接近于0, 模擬結(jié)果與理論一致． 注意到, 當(dāng)折線效應(yīng)不存在時(shí), 經(jīng)驗(yàn)P值接近名義上的顯著性水平0.05, 而隨著折線效應(yīng)增強(qiáng), 檢驗(yàn)功效增加, 并趨近于1, 這說(shuō)明本文提出的檢驗(yàn)統(tǒng)計(jì)量能夠有效識(shí)別折線效應(yīng)．

表4 K=2時(shí)各情形在不同折線效應(yīng)下統(tǒng)計(jì)量Tn均值及檢驗(yàn)功效的勢(shì)

本文亦使用所提方法“ESL．IND”“ESL．CHO”分析文獻(xiàn)[4]的縱向黃體酮數(shù)據(jù), 參數(shù)的估計(jì)值、 平均絕對(duì)誤差、 置信區(qū)間, 以及擬合曲線(結(jié)果見(jiàn)Github)． 實(shí)證結(jié)果顯示, 相較“SEG”“MKQR”方法, 所提估計(jì)具有明顯競(jìng)爭(zhēng)優(yōu)勢(shì)．

5 證明

為了證明定理1, 我們給出如下引理．

(16)

(16)式等價(jià)于

以及

(17)

(18)

結(jié)合(17)式及(18)式,

故而定理1可證．

定理2的證明記折點(diǎn)位置真值為τ0=(τ1, 0, …,τK,0)T． 由(9)式可得等式

經(jīng)計(jì)算,

定理3得證．

6 總結(jié)

本文基于指數(shù)平方損失提出縱向多折點(diǎn)回歸模型參數(shù)估計(jì)和統(tǒng)計(jì)推斷方法． 為處理折點(diǎn)回歸模型, 本文首先基于局部線性平滑方法將折點(diǎn)回歸模型轉(zhuǎn)化為普通的縱向線性模型． 然后為融合縱向數(shù)據(jù)中重復(fù)觀測(cè)間的相關(guān)性, 本文利用修正的Cholesky分解方法對(duì)重復(fù)觀測(cè)間的協(xié)方差陣進(jìn)行建模, 以提高回歸模型參數(shù)的估計(jì)效率． 本文討論了參數(shù)估計(jì)的大樣本性質(zhì), 并同時(shí)討論了指數(shù)平方損失函數(shù)中調(diào)諧參數(shù)的選擇、 折點(diǎn)個(gè)數(shù)的確定方法和折線效應(yīng)的檢驗(yàn)問(wèn)題等． 數(shù)值模擬和實(shí)證分析結(jié)果顯示本文所提方法可以為縱向多折點(diǎn)回歸模型的參數(shù)提供更為穩(wěn)健的、 有效的估計(jì)．