白永昕，錢曼玲，田茂再

（1.北京信息科技大學(xué) 理學(xué)院，北京 100192；2.墨爾本大學(xué) 數(shù)學(xué)與統(tǒng)計(jì)學(xué)院，澳大利亞 墨爾本3010；3.中國(guó)人民大學(xué)應(yīng)用統(tǒng)計(jì)科學(xué)研究中心，北京 100872；4.新疆財(cái)經(jīng)大學(xué) 統(tǒng)計(jì)與信息學(xué)院，烏魯木齊 830012；5.昌吉大學(xué) 數(shù)學(xué)與數(shù)據(jù)科學(xué)學(xué)院，湖南 昌吉 831100）

0 引言

隨著數(shù)據(jù)獲取技術(shù)的發(fā)展，在微陣列、蛋白質(zhì)組學(xué)、大腦圖像等領(lǐng)域都出現(xiàn)了超高維數(shù)據(jù)。在超高維數(shù)據(jù)中，協(xié)變量的維數(shù)可能遠(yuǎn)遠(yuǎn)大于樣本量，這給傳統(tǒng)的統(tǒng)計(jì)方法帶來(lái)了挑戰(zhàn)。高維數(shù)據(jù)通常是異質(zhì)的，協(xié)變量對(duì)條件分布中心的影響可能與他們對(duì)尾部的影響大不相同。因此，只關(guān)注條件均值函數(shù)可能會(huì)產(chǎn)生誤導(dǎo)。分位數(shù)回歸通過(guò)估計(jì)不同分位數(shù)水平上的條件分位數(shù)，能夠更完整地反映協(xié)變量和響應(yīng)變量之間的關(guān)系，而且分位數(shù)回歸對(duì)異常點(diǎn)有很強(qiáng)的魯棒性，在重尾分布下也會(huì)得到穩(wěn)健的估計(jì)。部分線性可加模型[1]作為一種半?yún)?shù)回歸模型，比參數(shù)模型更靈活，比非參數(shù)模型更有效。在復(fù)雜數(shù)據(jù)下研究部分線性可加分位數(shù)回歸模型的變量選擇問(wèn)題具有非常重要的理論和實(shí)際意義。

近年來(lái)，部分線性可加分位數(shù)回歸模型引起了學(xué)術(shù)界的廣泛關(guān)注。針對(duì)模型中線性部分的變量選擇問(wèn)題，一些研究者探索了不同的方法。例如，陳秀平和蔡光輝（2021）[2]使用非負(fù)Garrote方法選取重要變量；白玥和田茂再（2017）[3]、宋瑞琪等（2019）[4]系統(tǒng)對(duì)比了多種懲罰回歸方法（如Lasso、自適應(yīng)Lasso、SCAD、Elastic Net、組Lasso、組SCAD等）在不同自變量相關(guān)性和誤差項(xiàng)方差條件下的性能；Mazucheli等（2022）[5]則對(duì)線性分位數(shù)回歸模型的相關(guān)研究進(jìn)行了全面回顧。在高維數(shù)據(jù)分析場(chǎng)景中，盡管非凸懲罰，如SCAD 和MCP 具有更好的適應(yīng)性，但其復(fù)雜的性質(zhì)加大了優(yōu)化難度。為此，Wang 和Zhu（2016）[6]提出了適用于最小二乘回歸的arctan型懲罰，它相較于L0懲罰表現(xiàn)出更強(qiáng)的穩(wěn)定性及Oracle 性質(zhì)。Li 和Zhu（2008）[7]基于KKT 條件設(shè)計(jì)了Lasso 懲罰下的分位數(shù)回歸參數(shù)估計(jì)算法；而Wu和Lange（2008）[8]探討了中位回歸的快速貪婪坐標(biāo)下降算法；Wang等（2012）[9]進(jìn)一步將局部線性算法運(yùn)用到非凸懲罰的分位數(shù)回歸參數(shù)估計(jì)中，但該方法在高維協(xié)變量下的計(jì)算效率偏低。為應(yīng)對(duì)計(jì)算效率問(wèn)題，Peng 和Wang（2015）[10]研發(fā)了針對(duì)非凸懲罰的迭代坐標(biāo)下降算法（Iterative Coordinate Descent Algorithm，QICD），并證明了它的收斂性。模擬實(shí)驗(yàn)顯示，即使在極高維的情況下，QICD算法依然有效。在此基礎(chǔ)上，Sherwood 和Maidman（2022）[11]利用QICD算法對(duì)可加分位數(shù)回歸模型進(jìn)行了變量選擇和相關(guān)參數(shù)估計(jì)。

總體來(lái)看，現(xiàn)有文獻(xiàn)對(duì)部分線性可加分位數(shù)回歸模型變量選擇的研究已較為豐富，但仍存在一定的局限性：一是現(xiàn)有研究主要集中于協(xié)變量維度是固定的情況；二是并未同時(shí)考慮線性部分和可加部分的稀疏性。鑒于此，本文考慮了協(xié)變量維度發(fā)散情況下的部分線性可加分位數(shù)回歸模型，通過(guò)雙懲罰方法對(duì)模型中的線性部分和可加部分進(jìn)行變量選擇和穩(wěn)健估計(jì)，并推導(dǎo)了估計(jì)量的漸近性質(zhì)，以期推動(dòng)該問(wèn)題的研究進(jìn)展，獲得對(duì)部分線性可加分位數(shù)模型估計(jì)問(wèn)題更深入、更全面的理解。

1 部分線性可加分位數(shù)回歸模型

假定{(Yi，xi，zi):i=1，…，n} 是獨(dú)立同分布樣本，Yi是響應(yīng)變量。本文考慮協(xié)變量維度隨樣本量n變化的情況。記pn為隨樣本量變化的協(xié)變量維度，xi=(xi1，…，)和zi=(zi1，…，zid)分別是參數(shù)部分和非參數(shù)部分的協(xié)變量向量。給定(xi，zi),Yi的條件分位函數(shù)為=inf{t:F(t|xi，zi)≥τ}，其中，F(xiàn)(·|xi，zi)是Yi的條件函數(shù)。考慮如下部分線性可加分位數(shù)回歸模型：

2 部分線性可加分位數(shù)回歸模型的雙懲罰估計(jì)

2.1 基于Atan雙懲罰的估計(jì)

在實(shí)際數(shù)據(jù)中，通常并不知道哪些變量是重要變量。為了得到稀疏的估計(jì)量，最小化如下懲罰目標(biāo)函數(shù)：

其中，P(β，γ)表示關(guān)于參數(shù)β以及非參數(shù)函數(shù)gj的懲罰函數(shù)。為了構(gòu)造一個(gè)在保證模型稀疏性的同時(shí)還能保證估計(jì)函數(shù)光滑性的估計(jì)量，本文提出了Atan 雙懲罰函數(shù)：

2.2 漸近性質(zhì)

條件（1）是分位數(shù)回歸中使用的標(biāo)準(zhǔn)假設(shè)。條件（2）對(duì)于B樣條基函數(shù)是必要的，它可以用來(lái)有效地逼近滿足Holders條件的函數(shù)。條件（3）是一個(gè)可識(shí)別條件，對(duì)Oracle模型下的協(xié)變量和設(shè)計(jì)陣進(jìn)行了約束。假設(shè)xi四階矩有界就足夠了。條件（4）給出了用B樣條逼近非參數(shù)部分時(shí)樣條基的維度。條件（5）給出了β最小非零項(xiàng)個(gè)數(shù)的下限。條件（6）對(duì)Atan 雙懲罰中的調(diào)整參數(shù)λ1和α進(jìn)行了約束，類似的約束可參見(jiàn)文獻(xiàn)[6]。條件（7）限制了雙懲罰函數(shù)的一階導(dǎo)數(shù)和二階導(dǎo)數(shù)的變化速度。條件（8）用于證明Atan 估計(jì)的漸近正態(tài)性，與Lindeberg-Feller 中心極限定理中的Lindeberg條件有關(guān)。

定理1：假設(shè)正 則條件（1）至 條件（8）成立，則?η∈(0，1)，?常數(shù)C＞0，使得：

定理2：假設(shè)正則條件（1）至條件（8）成立，則：

3 算法

交替迭代算法的步驟如下：

在步驟2.2中，Atan懲罰很顯然是非凸函數(shù)，可以用它的一階泰勒近似代替非凸懲罰值，得到一個(gè)在當(dāng)前值β(t)處的凸目標(biāo)函數(shù)，即。本文使用分位迭代坐標(biāo)下降算法對(duì)步驟2.2中的目標(biāo)函數(shù)進(jìn)行最小化，限于篇幅，算法的具體步驟省略。

4 模擬研究

本文通過(guò)模擬研究對(duì)所提方法的性能進(jìn)行研究，并將其與現(xiàn)有方法（Lasso、SCAD和MCP懲罰）進(jìn)行比較。為了減輕計(jì)算負(fù)擔(dān)，將SCAD、MCP 以及Atan 雙懲罰中的參數(shù)分別設(shè)置為a1=3.7、a2=2.7 以及a3=0.005。對(duì)于樣本量和協(xié)變量維度，考慮n=200 和pn=100，500。在模擬研究中，重復(fù)模擬100次并考慮3個(gè)不同的分位點(diǎn)τ=0.3，0.5，0.7。

考慮如下模型：

對(duì)于參數(shù)估計(jì)的精度，通過(guò)以下指標(biāo)進(jìn)行評(píng)價(jià)：（1）均方誤差（MSE）：通過(guò)MSE 來(lái)衡量參數(shù)估計(jì)的精度。（2）在一個(gè)由均勻分布在[0,1]上的500 個(gè)點(diǎn)(t1，…，tT)組成的細(xì)網(wǎng)格上計(jì)算兩個(gè)分量函數(shù)的均方根誤差（RMSE），通過(guò)MSE來(lái)評(píng)估非參數(shù)函數(shù)的性能。在計(jì)算中，若參數(shù)估計(jì)值小于1e-06，則默認(rèn)其值為0。

對(duì)于隨機(jī)誤差項(xiàng)的分布，考慮如下三種不同的情形：（1）標(biāo)準(zhǔn)正態(tài)分布；（2）自由度為3的t 分布；（3）異方差的正態(tài)分布，即?i=Xi1ζi，其中，ζi～N(0，1) 且與Xi相互獨(dú)立。不同誤差分布下的模擬結(jié)果見(jiàn)表1至表3。

表1 隨機(jī)誤差項(xiàng)?～N(0，1)情況下的模擬結(jié)果

表2 隨機(jī)誤差項(xiàng)?～t(3)情況下的模擬結(jié)果

表3 隨機(jī)誤差項(xiàng)?～Xi1ζi， ζi～N(0，1) 情況下的模擬結(jié)果

從表1 至表3 可以看出，在不同隨機(jī)誤差項(xiàng)分布和不同pn值下，所有方法對(duì)非線性部分的擬合都比較相似。同時(shí)，在變量選擇上，這些方法的TPR都在1左右，說(shuō)明所有方法都可以選擇重要的變量。進(jìn)一步觀察發(fā)現(xiàn)，在大多數(shù)情況下，特別是當(dāng)隨機(jī)誤差項(xiàng)服從標(biāo)準(zhǔn)正態(tài)分布時(shí)，本文所提Atan 雙懲罰估計(jì)量的MSE 小于其他懲罰的估計(jì)量，這可能是因?yàn)锳tan 雙懲罰估計(jì)量是無(wú)偏的。同時(shí)，Atan 雙懲罰估計(jì)量不正確篩選重要變量的比例相比其他方法也更低。更重要的是，當(dāng)隨機(jī)誤差項(xiàng)服從t分布和異方差的正態(tài)分布時(shí)，其他懲罰估計(jì)量的TNR顯著下降，主要集中在0.5 左右，而Atan 雙懲罰估計(jì)量的TNR 保持在0.8左右。與此同時(shí)，其他懲罰估計(jì)量的FDR逐漸上升，而Atan 雙懲罰估計(jì)量的FDR 大部分保持在0.2 附近。不可忽視的是，本文所提方法以較大的MSE為代價(jià)，選擇了更精確的模型?？傊?，當(dāng)隨機(jī)誤差項(xiàng)服從重尾分布和異方差分布時(shí)，本文所提Atan雙懲罰的性能更好。

5 實(shí)證分析

將本文提出的Atan 雙懲罰方法應(yīng)用于一個(gè)包含315個(gè)癌癥篩查病人的血液樣本數(shù)據(jù)集。該數(shù)據(jù)集來(lái)源于http://lib.stat.cmu.edu/datasets/Plasma_Retinol，主要記錄了每個(gè)病人的β-胡蘿卜素血漿濃度、年齡、性別、體重、是否抽煙飲酒等14 個(gè)變量的數(shù)據(jù)。已有研究表明，血漿中的β-胡蘿卜素含量與患一些特定類型癌癥的風(fēng)險(xiǎn)相關(guān)，因此本文的目標(biāo)是找到影響β-胡蘿卜素血漿濃度的變量。Guo等（2013）[14]研究發(fā)現(xiàn)，Age（年齡）、Chol（每天攝入的膽固醇）和Fiber（每天攝入的動(dòng)植物纖維）與血漿β-胡蘿卜素水平存在非線性關(guān)系，其他變量則與血漿β-胡蘿卜素水平存在線性關(guān)系。因此，本文考慮采用部分線性可加分位數(shù)回歸模型對(duì)標(biāo)準(zhǔn)化后的數(shù)據(jù)集進(jìn)行變量選擇。表4 給出了不同懲罰下變量選擇的結(jié)果，其中，協(xié)變量分別為Quet（體重/身高的平方）、Calor（每天攝入的卡路里）、Fat（每天攝入的脂肪）、Alco（每周攝入的酒精）、Betad（每天攝入的β-胡蘿卜素）、Retd（每天攝入的視網(wǎng)醇）、Retpl（視網(wǎng)醇血漿濃度）、Sex（性別）、Smok（是否抽煙）、Vit（是否經(jīng)常使用維生素）。圖1給出了三個(gè)非參數(shù)成分在不同分位數(shù)下的估計(jì)曲線。從表4可以看出，在Lasso 懲罰下，協(xié)變量Fat和Retd未被選出，但是其他三個(gè)非凸懲罰均選出了這兩個(gè)協(xié)變量。同時(shí)，三種非凸懲罰的結(jié)果相似，但SCAD 和MCP 懲罰法高估了協(xié)變量Fat和Retd。Fat和Retd對(duì)變量選擇結(jié)果的影響相對(duì)較小。

圖1 不同分位數(shù)下年齡、攝入纖維、膽固醇的估計(jì)曲線

表4 不同懲罰下變量選擇的結(jié)果（τ=0.5）

為了進(jìn)一步評(píng)估本文方法的性能，將數(shù)據(jù)集隨機(jī)分為樣本量為215的訓(xùn)練集和樣本量為100的測(cè)試集。重復(fù)模擬100次并計(jì)算預(yù)測(cè)誤差。表5給出了100次模擬下不同篩選方法在0.3，0.5，0.7 分位點(diǎn)處選擇的平均變量數(shù)量以及預(yù)測(cè)誤差，反映了模型的平均復(fù)雜程度和預(yù)測(cè)能力；括號(hào)內(nèi)的值為相應(yīng)的標(biāo)準(zhǔn)誤，反映了不同方法在100次重復(fù)模擬中的波動(dòng)性。從表5 可以看出，在不同懲罰下，選出的模型往往不同。在Lasso懲罰下選出的模型比在非凸懲罰下選出的模型復(fù)雜程度更低，但同時(shí)預(yù)測(cè)精度也更低。對(duì)比三種非凸懲罰可以發(fā)現(xiàn)，本文提出的Atan 雙懲罰方法比SCAD 和MCP 懲罰預(yù)測(cè)精度更高，而且標(biāo)準(zhǔn)差也更小，即本文方法更穩(wěn)定。綜合來(lái)看，基于Atan雙懲罰的變量選擇結(jié)果較為理想，是一種相對(duì)穩(wěn)健的懲罰方法。

表5 模型擬合和預(yù)測(cè)

6 結(jié)論

本文研究了協(xié)變量維度pn發(fā)散且為超高維情況下的部分線性可加分位數(shù)回歸模型的變量選擇和穩(wěn)健估計(jì)問(wèn)題。首先，對(duì)于模型中的非參數(shù)函數(shù)，考慮用三次B 樣條函數(shù)進(jìn)行擬合。這種方法不僅在計(jì)算上十分便捷，而且通常能夠提供準(zhǔn)確的結(jié)果。為了實(shí)現(xiàn)超高維線性部分的稀疏性以及非參數(shù)函數(shù)的光滑性，本文提出一種Atan 雙懲罰估計(jì)量，并在一定的正則條件下推導(dǎo)了雙懲罰估計(jì)量的收斂速度和變量選擇的一致性。其次，為了解決所提方法中的優(yōu)化問(wèn)題，采用了一種迭代坐標(biāo)下降算法，即使在pn?n的情況下也能夠?qū)崿F(xiàn)快速收斂。模擬研究表明，當(dāng)隨機(jī)誤差項(xiàng)服從標(biāo)準(zhǔn)正態(tài)分布時(shí)，本文所提Atan 雙懲罰估計(jì)量的MSE 小于其他懲罰的估計(jì)量，因?yàn)锳tan 懲罰估計(jì)量是無(wú)偏的。值得注意的是，當(dāng)誤差分布存在重尾時(shí)，其他懲罰估計(jì)量的FDR 逐漸上升，而Atan 雙懲罰估計(jì)量的偽發(fā)現(xiàn)率仍能保持在較低水平，即Atan 方法選擇了更精確的模型。最后，將本文提出的方法應(yīng)用于一個(gè)包含315 個(gè)癌癥篩查病人的血液樣本數(shù)據(jù)的數(shù)據(jù)集。通過(guò)比較不同懲罰下的變量選擇結(jié)果發(fā)現(xiàn)，本文提出的Atan 雙懲罰方法在預(yù)測(cè)精度和標(biāo)準(zhǔn)誤方面均優(yōu)于SCAD 和MCP懲罰，表明該方法更穩(wěn)定?？傮w而言，基于Atan雙懲罰的選擇結(jié)果相對(duì)理想，是一種相對(duì)穩(wěn)健的懲罰方法。