基于自主學習與SCAD-Net正則化的回歸模型①

劉 杰, 陳浩杰

1(中國科學技術大學 管理學院, 合肥 230026)

2(中國科學技術大學 國際金融研究院, 合肥 230026)

基因組學研究的一個關鍵問題是如何確定與疾病相關的基因及其生物途徑, 常見的做法是通過將高維基因組數據(如微陣列基因表達數據)與各種臨床結果聯系起來構建疾病診斷預測模型. 然而, 迄今為止, 雖然許多基因生物標志物研究[1,2]已經完成, 但目前提出的眾多相關方法在臨床應用中均難以得到令人滿意的結果. 其原因主要在于研究樣本量太小[3,4], 從而導致統計效能降低, 進而得到可信度較低甚至錯誤的結論. 因此, 充足的樣本是產生有效統計分析和結論的必要條件. 另一方面, 數據收集技術的進步促使現行可用生物數據日益增多, 于是有學者提出了數據融合的思想, 即綜合多個數據集或有關結果. 然而, 盡管一些基因表達研究有著相同的目標, 但所用數據集通常是來自不同的處理設備、不同的數據平臺, 甚至彼此之間具有不同的數值尺度, 從而導致批量效應的存在. 因此, 直接整合不同的基因表達數據將會給統計分析帶來巨大挑戰.

為解決上述問題, 研究者們做了大量的工作, 主要分為以下兩類: 元分析和融合分析[5]方法. 元分析即利用統計的概念與方法去收集、整理以及分析之前學者針對某個主題所做的眾多實證研究. 然而, 元分析對一些必要條件較為敏感, 稍加違反就可能造成錯誤性結論[6]. 融合分析是對不同的數據集進行整合并以此作為研究數據集. 相比元分析, 融合分析具有更多的樣本從而更具統計效用. 近年來, 基于融合分析的方法層出不窮, 如Benito等[7]提出的距離加權判別法(DWD),Johnson等[8]提出的經驗貝葉斯方法(EB), Shabalin等[9]提出的跨平臺標準化方法(XPN), Deshwar等[10]提出的PLIDA方法以及Deng等[11]提出的WaveICA方法. 然而, 由于批量效應的存在, 且其來源復雜無法消除, 導致以上方法均可能給融合數據集帶來新的系統誤差, 使其變得更加復雜. 因此, 直接分析融合后的數據可能會引起一些問題[12,13], 需要提出一種新的方法來解決數據融合問題.

Kumar等[14]提出的自主學習(Self-Paced Learning,SPL)方法可以根據模型已經學習的內容自適應地識別簡單和困難樣本, 并且隨著模型訓練的不斷進行, 越來越多的困難樣本進入模型. SPL方法可以在很大程度上克服批量效應, 并且其應用較為廣泛, 目前已成功應用于各種機器學習問題[15]. 此外, Ma等[16]還對SPL方法的收斂性質進行了補充和討論, 使其在理論上更加豐富.

除樣本規模問題之外, 樣本維度是另一研究熱點.許多研究中的樣本維數通常遠遠大于樣本數量, 即常見的高維度低樣本問題. 這在生物統計中尤為常見, 如基因表達數據. 為解決該問題, 研究者們提出了許多正則化方法, 用于在回歸框架中識別與臨床表型相關的基因, 如Lasso[17]、SCAD[18]、Elastic-Net[19]、Fused Lasso[20]、Lars[21]、adaptive Lasso[22]、Group Lasso[23]以及L1/2+2混合正則化方法[24-26]. 然而, 以上正則化方法都存在共同的局限性, 即這些方法僅是從計算或算法的角度出發, 沒有利用任何先驗知識或信息.但對于許多復雜的疾病尤其是癌癥, 許多生物學途徑信息對于了解治療疾病具有較大的效用, 并且該信息可以從多年的生物醫學研究中獲得, 故將此種先驗信息納入模型考慮應該會有更好的預測效果.

本文將基于SPL方法構建一個更精確的基因表達預測模型. 首先我們將不同的基因表達數據集融合到一個統一的數據集中, 緊接著在線性回歸的背景下將SPL方法與SCAD網絡懲罰相結合得到最終的回歸預測模型. 具體來說, 該模型由3部分組成: (1) SCAD罰函數. 利用SCAD懲罰來增強模型的稀疏性, 該懲罰不僅為大系數提供了無偏估計, 并且具有較高的理論價值, 例如Oracle性質[18]； (2)基于網絡的懲罰, 利用網絡懲罰來實現基因調控網絡上相鄰節點系數之間的平滑； (3) SPL方法, 促使模型自適應地從簡單樣本(高置信度樣本)向復雜樣本(低置信度樣本)上過渡. SPL方法對于分析融合數據是至關重要的, 因為融合數據往往存在較大的噪聲以及異常值點.

本文接下來內容安排如下: 第1節提出了一個基于SCAD網絡懲罰的線性回歸模型, 緊接著介紹了自主學習(SPL)方法并將其與SCAD網絡懲罰相結合從而得到最終的預測模型； 第2節首先對SCAD網絡懲罰函數的理論性質進行簡單分析, 包括群組效應以及漸近性質； 然后給出一種求解本文所提出模型的有效算法； 在第3節中, 通過不同情形下的模擬數據以及在乳腺癌細胞系數據集上的分析結果來評估本文所提出模型的預測效果. 第4節是結論與展望.

1 SCAD網絡正則化與自主學習方法

1.1 SCAD網絡正則化

假設數據集D={(x1,y1),(x2,y2),···,(xn,yn)}, 其中xi=(xi1,xi2,···,xip)T表 示第i個 樣本,yi為對應的響應變量, 記X=(x1,x2,···,xp),Y=(y1,y2,···,yn)T. 進一步, 假設各個預測因子xi(i=1,2,···,p)經過標準化處理, 響應變量y經過去中心化處理, 從而有:

本文考慮最簡單的線性回歸模型:

式中, β=(β1,β2,···,βp)T為要估計的參數, εi表示均值為0, 方差為σ2的誤差項. 上述模型的平方損失函數可以表示為:

在許多研究當中, 樣本維度通常遠遠大于樣本數量,即高維度低樣本問題. 在這種情形下, 線性回歸并不能夠直接用來估計回歸參數. 由此, 引入了正則化方法, 即:

其中,P(β)表示正則化項. 高維變量選擇中常用的正則化方法為L1約束, 即Lasso方法, 具體可以表示為這里λ 表示任意非負數, 一般可使用k折交叉驗證方法確定. 由于L1罰函數具有奇異性,故基于L1懲罰的線性回歸模型可以將一些系數較小的參數壓縮為0從而達到變量選擇的效果. 但當 λ過大時, β估計量中系數較大的參數會存在較大偏差, 而當λ 過小時, β估計量則不夠稀疏. 為克服這一問題,Fan等[18]提出了SCAD懲罰函數, 其具體形式為:

其中,a為一個大于2的超參數, 根據文獻[18], 可將其設為3.7, 當然也可以通過交叉驗證的方法加以確定. 從表達式(1)可以發現, 當| β|較小時, 懲罰函數為線性函數；當| β|較 大時, 懲罰函數為二次懲罰； 當| β|很大時, 懲罰項為常數. SCAD懲罰函數關于 β的一階導函數為:

在非0處, 對任意 βj≈zj, 由二階泰勒展開可得:

此外, 對于協變量之間存在高相關性的問題, Zou等[19]提出了Elastic-Net懲罰函數, 其具體表達為Zeng等[27]提出了SCAD-L2懲 罰, 該懲罰同時結合SCAD和L2懲罰. 以上幾種方法都可以實現群組效應, 即具有強相關性的預測因子要么同時被選中, 要么同時被剔除. 然而, 上述幾種方法都是從計算的角度出發, 都沒有考慮先驗信息, 如網絡拓撲信息.

我們知道, 生物統計中的基因交互信息對于識別基因組模式具有重要價值. 該先驗信息可以用一個加權圖G=(N,E,W) 來表示, 其中N是網絡節點集, 表示p個預測因子,E={u～v}是圖中邊的集合, 表示節點u和v之間有邊相連,W={w(u,v)}則表示邊上的權重. 近年來, 網絡懲罰出現在大量現實應用中, 例如, Li等[28],Chen等[29]以及Wang等[30]利用基于網絡的L1懲罰對基因組數據進行回歸分析并進行變量選擇. 在這些研究當中, 網絡懲罰函數被定義為拉普拉斯矩陣的二次型. 然而, 在某些情況下,L1懲罰存在偏差并且可能導致結果不夠稀疏, 而SCAD懲罰可以避免過度懲罰并且具有良好的統計性質. 因此, 本文給出基于網絡的SCAD-Net懲罰函數(SCAD Network-based penalized fuction, SCAD-Net):

其中,L表示拉普拉斯矩陣, 根據文獻[31], 可將其定義為:

其中, 參數 λ1與 λ2分別控制參數估計的稀疏度與光滑度. 進一步, 在線性回歸背景下, 可以得到基于SCADNet懲罰的線性回歸模型(SCAD-Net penalized Linear regression, SNL):

其中, 第1項表示線性回歸的損失函數； 第2項表示SCAD懲罰函數, 保證參數估計的稀疏性, 并且強化結果的可解釋性； 最后一項表示基于網絡的懲罰函數, 保證參數估計的光滑性, 并且將網絡結構信息與基因表達數據相融合.

1.2 自主學習策略

盡管正則化方法在基因數據分析, 變量選擇等方面扮演非常重要的角色, 但最終得到的結論卻鮮少在臨床中得到應用. 這是因為上述結論都是基于小樣本數據, 導致結果的可信度較低. 為解決這一問題, 有人提出通過整合不同的數據集來生成人工大樣本數據.然而, 這些數據整合的方法并不能消除內部偏差, 甚至可能會增加新的誤差.

受人類學習機制的啟發, Kumar等[14]提出了自主學習(SPL)方法, 該方法首先通過低噪聲樣本學習一個基礎模型, 然后通過高噪聲樣本學習使模型變得更加穩健, 該方法可以顯著提高融合數據集的統計分析效用. 并且Kumar表示, 通過引入一個懲罰項, 可以將自主學習方法視為優化模型, 具體可以表示為:

為了加強對融合數據分析的準確性與魯棒性, 本文將SPL方法與SCAD-Net正則化在線性回歸的背景下相結合, 從而得到最終的回歸模型(Self-paced learning and SCAD-Net penalized Linear regression, SSNL):

其中, 第1項表示加權的線性回歸模型, 最后兩項表示SCAD-Net懲罰函數.

2 理論性質及求解方法

2.1 理論性質

本小節我們給出與SCAD-Net正則化方法相關的性質, 包括群組效應以及在p固 定且n→∞ 情形下的漸近性質.

2.1.1 群組效應

Huang等[32]證明了SCAD-Net懲罰函數具有群組效應, 如引理1與引理2所示, 其具體證明過程見文獻[32].

引理1可看作Zou等[19]中引理2的進一步結果, 其保證在兩個預測變量相等時, 估計參數具有群組效應.

則有:

引理2給出SCAD-Net懲罰函數群組效應的量化描述, 即在滿足以上條件的前提下, 兩個參數的差異具有上界約束. 進一步, 若樣本相關系數 ρ趨于1, 則兩個估計參數幾乎相同.

2.1.2 漸近性

SCAD-Net懲罰線性回歸的目標函數為:

證明: 定義

進一步, 根據式(2)和式(3)可知:

對任意 βj≈zj, 有:

同樣地, 關于第三項有:

因此, 在有限維收斂的情況下有:

又Vn為凸函數且V有最小值, 可得:

證畢.

2.2 求解方法

本節給出模型SSNL的求解算法, 具體如下:

(1) 固定v更新 β時, 相當于解決SCAD-Net正則化回歸問題, 即SNL. 本文利用坐標下降法進行求解,具體來說, 式(8)關于 βj(j=1,2,···,p)求導, 可得:

令其等于0, 有:

具體更新算法如算法1所示.

算法1. SNL t=0,βj(t)=βj0 β0 1. 令 , 其中 為OLS估計量, 各調優參數可利用交叉驗證方法得到；βj(t),j=1,2,···,p 2. 通過式(9)依次更新 ；t←t+1 3. 令 , 重復步驟2直至:p∑j=1|βj(t)-βj(t-1)|＜10-6

(2) 固定 β更新v時, 式(8)關于vi求導, 可得:

進而有:

對于樣本i, 若其損失小于超參數 τ, 則可將其視為高質量樣本, 相對應的vi設為1, 否則設為0. 顯然, 對于樣本損失小于 τ的樣本會被納入模型中. 一旦得到v,我們進一步增大 τ的值, 這樣具有更大損失的樣本將會進入模型當中, 重復上述步驟直至收斂, 完整算法如算法2所示.

3 數值結果

3.1 模擬數據分析

為檢驗本文所提出SSNL模型的預測表現, 我們首先按照以下方式模擬出一個簡單的基因調控網絡:假設有200個轉錄因子(TFs), 每個轉錄因子調控10個基因,由此產生由2200個基因(節點)組成的生物基因調控網絡, 轉錄因子之間以及與其調控的基因之間形成網絡的邊. 為了簡單起見, 我們進一步假設模型中只有4個轉錄因子以及其調控的基因與響應變量y有關. 對于第一個模型, 我們按照以下方式來生成相關數據:

(3) ε為誤差項, 且εi～5×N(0,1).

(4) 200個轉錄因子服從標準正態分布, 即xTFj～N(0,1),j=1,2,···,200.

(5) 每個TF與其調控的單個基因均服從二元正態分布, 且相關系數為 ρ.

對于模型2, 我們假設

其他設置與模型1完全一致. 該模型假設同一轉錄因子所調控的基因既可以對響應變量y產生正影響,也可以對其產生負影響.

我們將模擬數據分為訓練集和測試集, 其中訓練集所占比例為70%. 在實際應用中, 人們通常使用k(k=3,5或10)折交叉驗證的方法來選擇調優參數, 然而, 不同的k折交叉驗證的方法很可能產生非常相似的預測結果[33,34]. 此外, 研究結果表明將交叉驗證的折數從10減少到3, 可以使算法的計算時間減少一半以上. 鑒于本文所提出的模型存在多個調優參數, 因此, 最終我們采用3折交叉驗證的方法來選擇最佳調優參數.

將基因相關系數 ρ分別設為0.2、0.5以及0.7. 每種情況都獨立重復模擬50次, 并計算得到相應的均方誤差(PMSE). 我們還進一步給出表征模型特征選擇效果的兩個指標, 分別是P和TP. 其中P表示模型中非零系數的個數, TP表示實際模型中非零系數的個數.表1給出了各模型在不同情形下的模擬結果, 其中Lasso-Net表示Lasso和網絡正則化； SCAD-Net表示不使用自主學習方法的SCAD和網絡正則化.

如表1所示, 在不同模型不同相關系數下, 本文提出的SSNL方法相比于Lasso-Net和SCAD-Net方法均給出最小的PMSE. 此外, 在識別真正相關基因時,SSNL相比于其他兩種方法具有最高的準確性(TP). 例如, 在 ρ=0.7 的情形下, 利用SSNL模型計算得到的TP值均超過40, 幾乎達到模型的真實值44. 上述結果表明SSNL方法在處理高維度低樣本、高噪聲、高相關性的復雜數據集時具有良好的表現.

表1 各模型在不同情形下的模擬結果

3.2 實際數據分析

為進一步論證SSNL模型的預測效果, 我們收集得到了乳腺癌細胞系數據集. 該數據集共有56個樣本,其中每個樣本都隸屬于一個確定的細胞亞型, 通過對其進行編碼可以得到一個數值型響應變量. 此外, 每個乳腺細胞樣本包含39 653個基因, 并且這些基因之間存在交互關系. 通過加權基因共表達網絡分析, 我們可以得到相應的加權網絡. 然后將基因表達數據與該調控網絡相結合, 得到最終的研究數據集. 我們旨在探索基因網絡與關注的表型之間的關聯關系以及網絡中的核心基因.

我們將數據集隨機打亂, 使約70%的樣本成為訓練樣本, 剩余30%的樣本作為測試樣本. 類似于上文模擬中的情形, 我們采用3折交叉驗證來估計得到最佳的調優參數. λ1與 λ2的候選值均來自于{0.01: 0.1:5}(起始值: 步長: 終值), μ來自于{1.1: 0.1: 3}以及 τ來自于{0.1: 0.05: 0.5}. 獨立重復10次, 計算得到相應的均方誤差(PMSE)以及模型中非零系數的個數P, 具體結果如表2所示.

表2 各模型在乳腺癌細胞系數據集上的結果

從表2可以看出, 本文提出的SSNL方法給出了最小的PMSE, 其表現顯著優于Lasso-Net方法, 且優于不使用自主學習的SCAD-Net方法. 此外, 在特征選擇方面, 盡管3種方法的數值表現效果相當, 但SSNL方法仍優于其他兩種對比方法. 上述結果再次說明本文所提出的SSNL模型在處理高維復雜網絡數據集時具有良好的表現.

4 結論與展望

融合分析為基因組研究提供了一種有效的分析角度. 傳統的融合分析方法是將多個數據集組合成一個集成的數據集, 然后直接對數據進行分析. 然而, 這種集成方法非但不能消除內部偏差, 甚至可能給融合數據集增加新的隨機噪聲和估計誤差, 從而降低融合分析的統計功效. 本文提出了一種新的融合分析模型SSNL, 該模型融合了自主學習(SPL)和SCAD-Net正則化方法. 一方面, SPL方法能夠先從低噪聲樣本中學習出一個基本模型, 然后通過高噪聲樣本學習使得模型更加穩健. 另一方面, 特征選擇是SSNL模型的重要組成部分. SCAD罰函數是一種常見的特征選擇方法,但SCAD罰函數僅是從計算的角度出發, 沒有利用任何先驗信息. 故在已有研究的基礎上, 本文給出了結合網絡結構信息的SCAD-Net懲罰, 并對這一問題進行了一些理論探究, 包括群組效應和漸近性質. 不同情形下的模擬分析結果以及在乳腺癌細胞系數據集上的分析結果均表明, SSNL方法在處理高維復雜網絡數據集時具有良好的預測表現.

本文使用3折交叉驗證(CV)方法來選擇SSNL模型中出現的懲罰參數. 然而, 當遇到多個超參數時,使用CV方法進行網格搜索需要消耗大量的時間與內存. 最近, 一種進化計算(EC)方法被用來調整懲罰參數, 并且表現良好[35]. 針對本文情形, EC方法可能是一個更好的選擇. 此外, 我們還考慮將SPL+SCADNet方法拓展到其他回歸模型中, 如廣義線性回歸等.