次序L1 懲罰估計在高維數(shù)據(jù)分析中的應用

郭依菏，訾雪旻

（天津職業(yè)技術師范大學理學院，天津300222）

隨著生物信息學、生物醫(yī)學、計量經(jīng)濟學和機器學習等領域的迅猛發(fā)展，高維數(shù)據(jù)在實際應用中越來越受到關注。為了能夠全面地分析數(shù)據(jù)，研究者總是盡可能地收集與研究問題相關的變量。如醫(yī)學研究人員在研究哪些基因變異會導致某種病變時，需要同時追蹤成千上萬的基因序列的表達。但是基于檢測成本等客觀因素，研究人員往往只能得到很小的樣本量。對于變量的數(shù)量遠遠大于觀測值的數(shù)量這種情況，尋找合適的變量選擇方法，提高變量篩選的效率及正確率，成為近年來統(tǒng)計學界的熱點問題之一。Foster 等[1]提出風險膨脹準則（risk inflation criterion，RIC），使用L0范數(shù)作為懲罰項，但對于高維情況，此方法通常會導致模型中包含很多不相關的變量，特別是在稀疏情況下，預測能力表現(xiàn)相當差。Tibshirani[2]提出套索算法（least absolute shrinkage and selection operator，LASSO），其通過使用L1范數(shù)作為懲罰項，克服了傳統(tǒng)方法在變量選擇上的不足，成為較流行的處理高維模型中變量選擇的方法。Bogdan 等[3]受 Benjamini-Hochberg 多重檢驗校正方法的啟發(fā)，在LASSO 方法的基礎上，提出了一種新的變量選擇方法——次序L1懲罰估計（Sorted L1Penalized Estimation，SLOPE）方法[4]。本文利用數(shù)值模擬方法分析比較了SLOPE 和LASSO 方法，結果表明SLOPE 方法對一般稀疏回歸模型有很好的篩選效果，同時利用實證分析研究了SLOPE 方法的有效性。

1 次序L1 懲罰估計算法

設n 維響應向量y 是由一個簡單的線性模型生成的

式中：Xn×p為設計陣；β 為 p 維向量回歸系數(shù)；Z 為 n×1的隨機誤差項。

首先考慮在正交設計下，令X′X=Ip，且隨機誤差項Z 獨立同分布于正態(tài)分布即N（0，σ2In）。此時不需考慮交互效應，定義?為響應變量y 和設計陣X 的乘積，則回歸模型可以改寫為

為了控制錯誤發(fā)現(xiàn)率水平在q∈[0，1]，首先采用BH方法對響應變量以降序排序，即令與之對應的假設檢驗分別為 H(1)，…，H(p)，定義參數(shù)iBH如下

當i≤iBH時，拒絕原假設H(i)。集合i 為空時，定義iBH=0。與傳統(tǒng)的Bonferroni 校正方法相比，BH 方法是一種自適應過程，因為拒絕的閾值iBH由數(shù)據(jù)直接驅(qū)動。SLOPE 方法選取λBH（i）=Φ-1（1-q）作為閾值，即標準正態(tài)分布第（i·q/2p）百分位數(shù)，其中 qi= i·q/2p。由于打破了傳統(tǒng)的單一懲罰，通過BH 方法確定的自適應閾值，對真實信號的稀疏性具有敏感性及適應性，且具有漸進最優(yōu)性質(zhì)。

類似于LASSO 方法，SLOPE 方法的解可以轉(zhuǎn)化為如下優(yōu)化問題

由式（4）得到的b 項是由相應的λi確定的，滿足|b|(1)≥|b|(2)≥…≥|b|(p)，b 項相當于次序 L1范數(shù)。文獻[10-12]中討論了次序L1范數(shù)的凸性，相關的應用參見文獻[13-15]。正如次序L1范數(shù)相當于L1范數(shù)的擴展，SLOPE 方法可以看作LASSO 方法的擴展。不同之處在于SLOPE 方法中不同的變量根據(jù)其相對重要性受到不同程度的懲罰，即對值較大的系數(shù)的懲罰比對值較小的系數(shù)的懲罰更大。在式（4）中，雖然SLOPE 方法在建立過程中借鑒了BH 方法的思想，但SLOPE 方法不是使用一個標量閾值來求解，因此即使在正交變量的情況下，這些過程也不是完全等價的。

2 在一般稀疏回歸模型中的應用

2.1 在多重檢驗中的應用

本節(jié)展示了在一個具有特定相關結構的多重檢驗問題中，如何將SLOPE 作為一個有效的多重比較控制程序?？紤]以下情況，醫(yī)學家在5 個隨機選擇的實驗室中各進行了1 000 個實驗，得出的觀察結果可以建模為

其中實驗室處理效應τj獨立同分布于正態(tài)分布隨機誤差項 zi，j獨立同分布于正態(tài)分布且 τ 和 z 序列之間相互獨立。相應的假設檢驗原假設為 H0，j：μj=0，則雙邊備擇假設為 H1，j：μj≠0。對所有 5 個實驗室的分數(shù)進行平均，結果是

其中 ε～N（0，Ip），問題轉(zhuǎn)換為經(jīng)典的變量選擇問題：識別回歸系數(shù)向量μ 的非零組成部分，注意到Σ不是對角對陣，而Σ-1/2是對角占優(yōu)的。如當σ2=1，ρ=0.5 時，解得通過選擇Σ-1/2作為設計陣，可以將強正相關多重檢驗問題近似于正交設計下變量選擇的問題，此時可以運用SLOPE 進行變量選擇。

2.2 λ的選取

當設計陣的列是隨機的相互獨立的觀測數(shù)據(jù)時，給變量選擇的回歸模型中引入了一系列復雜的因素，由于其具有隨機性，列之間的內(nèi)積不等于0。而偏離正交的情況會影響L1懲罰估計結果的一致性。SLOPE通過BH 方法確定的自適應閾值λ 是遞減序列，類似于一個逐步下降過程，可以通過逐步修改合并λ 序列，削減由于相關變量收縮而產(chǎn)生的額外干擾。

類似于LASSO 方法中λ 序列的選擇，假設X 的列有一個單位范數(shù)，且Z～N（0，1）。首先給出了LASSO方法的最優(yōu)性條件

由 y=Xβ +Z 得

式中：ηλ為軟閾值算子，ηλ（t）=sgn（t）（|t|-λ）+。記vi=則

記 S 為 βi的支撐集，XS為滿足 βi≠0 的相關變量的子集，βS表示相應系數(shù)的值。根據(jù)式（10），為了找到Xi與非零集合之間的相關關系，考慮表達式假設βS和λ 的值的大小滿足對任意正整數(shù)在不失一般性的條件下，進一步假設βj≥0。由于LASSO 方法中滿足 KKT 條件

而對于 SLOPE 方法，λ 是一個非增序列（λ1，…，λp），而不是一個定值。將λ 序列代入上式，非零部分的估計值為

式中：λS=（λ1，…，λk）′，故由式（13）可得

在高斯設計下，設X 獨立同分布于標準正態(tài)分布N（0，1/n），i?S，則

在 SLOPE 方法中，λ 是一個遞減序列，其中 λBH（1）可以簡單看作用于選取第一個進入模型的變量，選取λBH（1）=λG（1），而在下一步，由于選取λG（1）造成了回歸系數(shù)的輕微收縮，需要考慮方差的輕微增加，故使用λG（2）代替λBH（2）

類似的，接下來選取λ3時需要考慮選取λG（2）和λG（2）造成回歸系數(shù)的收縮，得到遞推公式為

3 SLOPE 方法與LASSO 方法的比較

LASSO 方法通過使用L1范數(shù)作懲罰項，其本質(zhì)是在最小二乘估計中簡單地應用軟閾值規(guī)則，選擇的是一個固定的標量閾值。這相當于將所有p 值與一個固定的閾值λ 進行比較，轉(zhuǎn)化為如下優(yōu)化問題

利用L1正則化項在零點的奇異性，將一部分不重要的變量系數(shù)的回歸系數(shù)壓縮到0。但由于使用L1范數(shù)造成的收縮，使得預測效果存在一定的偏差且預測精確度較低。

SLOPE 方法用排序后的L1范數(shù)作為懲罰項，相當于一個凸優(yōu)化程序，在對回歸參數(shù)進行估計的同時也實現(xiàn)了對變量的篩選，而且對參數(shù)估計比較大的變量回歸系數(shù)也趨于增大。 但SLOPE 方法的思想與LASSO 完全不同，SLOPE 方法通過BH 方法得到的自適應的閾值對數(shù)據(jù)的稀疏性和大小很敏感，在控制FDR 達到名義水平的同時可以適應未知數(shù)據(jù)的稀疏性。同時通過修正λ 序列來平衡偏離正交性帶來的影響，并可以運用到一般稀疏回歸模型中進行變量選擇。

接下來利用數(shù)值模擬方法分析比較SLOPE 和LASSO 方法擬合稀疏回歸模型。首先，生成 X1000×1000為稀疏設計陣，其中設置非零變量數(shù)k 分別為10，20，…，50，令FDR 名義水平為q=0.1 時，運用R 程序分別使用SLOPE 和LASSO 方法進行變量選擇，得到Lasso 和SLOPE 比較結果如表1 所示。

表1 LASSO 和SLOPE 比較結果

由表1 可以看出，SLOPE 和LASSO 方法都可以控制FDR 在名義水平下。其中對于非零特征數(shù)較少的情況，SLOPE 方法可以把所有非零元素準確識別出來，在非零元素個數(shù)增加時，在近百次的實驗中分別出現(xiàn)了少數(shù)誤判的情況，此時TPR 為0.98，預測效果較好，有較高的準確率。而LASSO 方法只有在非零元素較少時表現(xiàn)良好，隨著非零元素的增多，無法適應數(shù)據(jù)的稀疏性，如非零元素為50 的情況，其平均只識別出14 個元素，控制FDR 在名義水平下，但TPR 只有0.28，無法在錯誤發(fā)現(xiàn)率和正確發(fā)現(xiàn)率中達到平衡，預測效果存在一定的偏差。相比之下，在一般情況下SLOPE 方法擬合稀疏回歸模型效果更好，是一種準確且有效的變量選擇方法。

4 實證分析

本文從Broad Institute 癌癥項目數(shù)據(jù)集中選取了149 名乳腺癌患者，篩選出1 213 種與乳腺癌相關的基因表達數(shù)據(jù)作為樣本①。記n=149，1 213 種基因分別記為 V1，V2，…，V1213，此時變量的數(shù)量遠遠大于觀測值的數(shù)量，需要采取變量選擇方法進行降維。分別運用SLOPE 和LASSO 方法來篩選數(shù)據(jù)，將一些無意義或者意義極小自變量的回歸系數(shù)壓縮至0，進行變量選擇。設置FDR 名義水平為q=0.1，通過SLOPE 和LASSO 方法分別得到了46 和19 個非零回歸系數(shù)，相應篩選出SLOPE 和LASSO 方法得到的顯著基因變量如表2 所示。

為驗證這2 種變量選擇方法的有效性，通過主成分分析對數(shù)據(jù)集進行降維處理，得到的碎石圖如圖1所示。從圖1 可以看出，折線圖在43 個成分之后趨于平緩，即至少需要43 個變量才能包含90%以上的信息。

表2 SLOPE 和LASSO 方法得到的顯著基因變量

圖1 主成分分析得到的碎石圖

由表2 可知，LASSO 方法篩選出19 個顯著變量，SLOPE 方法篩選出46 個顯著變量，通過這2 種方法選擇出來的相同變量有15 個。結合主成分分析結果，選取19 個變量最多只能包含原始變量65%左右的信息，選取46 個變量可以保留原始變量95%以上的信息。因此，可以推出相較于LASSO 方法過度壓縮數(shù)據(jù)，SLOPE 方法可以更好地適應未知數(shù)據(jù)的稀疏性，有較高的預測準確性?？刂苢=0.1 時，SLOPE 方法得到的回歸系數(shù)估計結果如表3 所示，其中對于相對更為顯著的變量加上更大的回歸系數(shù)。

表3 SLOPE 方法得到的回歸系數(shù)估計結果

表3 中回歸系數(shù)的均值μ=-15.76，標準差σ=99.84，因此得到的回歸系數(shù)具有顯著性差異。且在SLOPE 方法選擇出的46 個變量中，與LASSO 方法選擇出的相同變量有15 個，此15 個變量用SLOPE 方法得到的相應的回歸系數(shù)均為較大的值。由于SLOPE 方法對于更為顯著的變量回歸系數(shù)更大這一特性，可知在控制FDR 在名義水平q=0.1 時，LASSO 方法過于保守，只能篩選出較為顯著的一部分變量，由于使用的是一個固定的閾值，無法很好地適應未知數(shù)據(jù)的稀疏性，導致了過度壓縮數(shù)據(jù)。而SLOPE 方法使用的是基于BH 方法選擇的自適應的序列，可以適應數(shù)據(jù)的未知稀疏性，在控制整體的錯誤發(fā)現(xiàn)率的同時，達到了增加勢、控制第一類錯誤概率的目的，有更好的篩選效果。

5 結 語

本文基于SLOPE 方法研究對于高維數(shù)據(jù)的變量選擇效果，通過數(shù)值模擬表明，相對于傳統(tǒng)的變量選擇方法，SLOPE 可以在控制整體的錯誤發(fā)現(xiàn)率的同時，適應未知數(shù)據(jù)的稀疏性。但使用SLOPE 分析超高維數(shù)據(jù)時仍存在算法運行緩慢、結果不夠精確等問題，今后的工作可以結合交叉驗證或樣本拆分等方法提高方法的準確性和運行效率。

注 釋：

① 數(shù)據(jù)來源 http：//portals.broadinstitute.org/cgi-bin/cancer/datasets.cgi。