超高維數(shù)據(jù)下部分線性可加分位數(shù)回歸模型的變量選擇

白永昕，錢曼玲，田茂再

（1.北京信息科技大學(xué) 理學(xué)院，北京 100192；2.墨爾本大學(xué) 數(shù)學(xué)與統(tǒng)計學(xué)院，澳大利亞 墨爾本3010；3.中國人民大學(xué)應(yīng)用統(tǒng)計科學(xué)研究中心，北京 100872；4.新疆財經(jīng)大學(xué) 統(tǒng)計與信息學(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ù)可能遠遠大于樣本量，這給傳統(tǒng)的統(tǒng)計方法帶來了挑戰(zhàn)。高維數(shù)據(jù)通常是異質(zhì)的，協(xié)變量對條件分布中心的影響可能與他們對尾部的影響大不相同。因此，只關(guān)注條件均值函數(shù)可能會產(chǎn)生誤導(dǎo)。分位數(shù)回歸通過估計不同分位數(shù)水平上的條件分位數(shù)，能夠更完整地反映協(xié)變量和響應(yīng)變量之間的關(guān)系，而且分位數(shù)回歸對異常點有很強的魯棒性，在重尾分布下也會得到穩(wěn)健的估計。部分線性可加模型[1]作為一種半?yún)?shù)回歸模型，比參數(shù)模型更靈活，比非參數(shù)模型更有效。在復(fù)雜數(shù)據(jù)下研究部分線性可加分位數(shù)回歸模型的變量選擇問題具有非常重要的理論和實際意義。

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

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

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

假定{(Yi，xi，zi):i=1，…，n} 是獨立同分布樣本，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ù)回歸模型的雙懲罰估計

2.1 基于Atan雙懲罰的估計

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

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

2.2 漸近性質(zhì)

條件（1）是分位數(shù)回歸中使用的標(biāo)準(zhǔn)假設(shè)。條件（2）對于B樣條基函數(shù)是必要的，它可以用來有效地逼近滿足Holders條件的函數(shù)。條件（3）是一個可識別條件，對Oracle模型下的協(xié)變量和設(shè)計陣進行了約束。假設(shè)xi四階矩有界就足夠了。條件（4）給出了用B樣條逼近非參數(shù)部分時樣條基的維度。條件（5）給出了β最小非零項個數(shù)的下限。條件（6）對Atan 雙懲罰中的調(diào)整參數(shù)λ1和α進行了約束，類似的約束可參見文獻[6]。條件（7）限制了雙懲罰函數(shù)的一階導(dǎo)數(shù)和二階導(dǎo)數(shù)的變化速度。條件（8）用于證明Atan 估計的漸近正態(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ù)，可以用它的一階泰勒近似代替非凸懲罰值，得到一個在當(dāng)前值β(t)處的凸目標(biāo)函數(shù)，即。本文使用分位迭代坐標(biāo)下降算法對步驟2.2中的目標(biāo)函數(shù)進行最小化，限于篇幅，算法的具體步驟省略。

4 模擬研究

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

考慮如下模型：

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

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

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

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

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

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

5 實證分析

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

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

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

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

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

6 結(jié)論

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