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

白永昕，錢曼玲，田茂再

（1.北京信息科技大學 理學院，北京 100192；2.墨爾本大學 數學與統計學院，澳大利亞 墨爾本3010；3.中國人民大學應用統計科學研究中心，北京 100872；4.新疆財經大學 統計與信息學院，烏魯木齊 830012；5.昌吉大學 數學與數據科學學院，湖南 昌吉 831100）

0 引言

隨著數據獲取技術的發展，在微陣列、蛋白質組學、大腦圖像等領域都出現了超高維數據。在超高維數據中，協變量的維數可能遠遠大于樣本量，這給傳統的統計方法帶來了挑戰。高維數據通常是異質的，協變量對條件分布中心的影響可能與他們對尾部的影響大不相同。因此，只關注條件均值函數可能會產生誤導。分位數回歸通過估計不同分位數水平上的條件分位數，能夠更完整地反映協變量和響應變量之間的關系，而且分位數回歸對異常點有很強的魯棒性，在重尾分布下也會得到穩健的估計。部分線性可加模型[1]作為一種半參數回歸模型，比參數模型更靈活，比非參數模型更有效。在復雜數據下研究部分線性可加分位數回歸模型的變量選擇問題具有非常重要的理論和實際意義。

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

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

1 部分線性可加分位數回歸模型

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

2 部分線性可加分位數回歸模型的雙懲罰估計

2.1 基于Atan雙懲罰的估計

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

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

2.2 漸近性質

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

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

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

3 算法

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

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

4 模擬研究

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

考慮如下模型：

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

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

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

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

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

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

5 實證分析

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

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

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

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

表5 模型擬合和預測

6 結論

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