基于預測變量圖結構的高維邏輯回歸模型*

黃文靜， 鄧 丹， 杜杰琳， 吳明月

(重慶醫科大學 公共衛生與管理學院，重慶 400016)

0 引 言

隨著信息和技術的快速發展，數據搜集能力越來越強，數據的規模和復雜程度都是前所未見的，特別是激增的變量維數給傳統的統計學帶來了巨大挑戰。在高維情形下，變量個數p要大于樣本量n，而且高維變量之間通常具有很強的相依性，此時，設計矩陣X不是滿秩的，回歸參數的估計往往也不是唯一的[1]。因此，如何對高維變量建模和參數估計是統計學關注的重點問題。正則化方法是處理高維數據常用的一種方法。Lasso[2]是人們熟知的正則化方法，但是當數據中存在某種結構(如組結構)時，Lasso不能很好地利用這種結構。因此，作為Lasso的推廣，隨后進一步提出Group Lasso[3]和Overlapping Group Lasso[4]等，這些方法在高維數據分析中均得到了廣泛的應用，但是這些方法通常需要提前給出預測變量的組結構，然而在有些應用中很難得到這種結構。相比組結構，預測變量的圖結構更容易獲得。在生物學研究中，每個個體包含了成千上萬的基因，這些基因有的單獨作用，有的相互作用，而基因之間相互作用的大量信息可以用來構建預測變量的圖結構，其中基因表示圖中的節點，調控關系表示圖中的邊。即使在沒有先驗信息的條件下，仍然可以通過協差陣(精度矩陣)的稀疏估計[5-7]構造這些基因的圖結構。Yu和Liu[8]在線性模型的框架下,提出一種利用預測變量圖結構的方法，該方法具有一定的普遍性。目前，邏輯回歸是處理分類數據的有效工具，正則化方法也被廣泛應用于處理大p小n的邏輯回歸問題[9-10]，本文將在文獻[8]的基礎上將圖結構應用到邏輯回歸模型中，為了評價這一方法的性能，將其與現有的許多方法進行比較，研究比較了不同圖結構下的仿真實例，還將其應用到乳腺癌基因實例中，仿真實驗結果和實際數據分析表明：該方法在估計、預測和模型的變量選擇上均具有優越性。

1 基于預測變量圖結構的邏輯回歸模型

1.1 邏輯回歸模型的改進

設(X,Y)為隨機變量，其中X∈Rp，Y∈{0,1}。假設X服從均值為Op×1，協差陣為Σ的多元正態分布。給定X=x，Y的條件分布記為Pβ*(Y=1|x)=pβ*(x)。邏輯回歸模型為

損失函數為

(

β

)=

(

β

;

x

,

Y

)= -

Yf

β

(

x

)+ log(1+exp(

f

β

(

x

)))

相應的風險函數和經驗風險函數分別記為

P(β0,β)=E(β0,β;X,Y)

Pn(β0,β)=

由文獻[12]的定理2.1及文獻[13]的條件3.1，有

β*∝Σ-1(η(y)-μ)

(1)

令Σ-1=Ω=(ωij)i,j=1,2,…,p，其中Ω為精度矩陣，由式(1)可知，β*可表示為β*=Ωγ，γ=η(y)-μ。因此，文獻[10]的方法可以推廣到邏輯回歸模型，于是有

(2)

(3)

基于預測變量的圖結構G確定領域N1,N2,…,Np(注意到j∈Nj,j∈Nj,j=1,…,p)。

求解如下優化問題：

(4)

1.2 優化求解及權重參數的選取

(5)

目前，有很多R包可以求解優化問題式(5)，如grpLasso[9]，grpreg[15]以及gglasso[16]。本文利用Zeng和Breheny[9]提供的R包求解優化問題式(5)，該方法能夠直接求解重疊組結構的Group Lasso 問題。

相比于目前流行的正則化方法(如Lasso，Group Lasso等)，本文的方法更具一般性。當預測變量的圖結構沒有邊時，本文方法等價于Adaptive Lasso；當預測變量的圖結構由若干個獨立的完全圖組成，本文方法等價于Group Lasso；當預測變量的圖結構為完全圖時，本文方法與嶺回歸有相同的非零解[11]。因此，本文的方法是一種更加一般的方法，Adaptive Lasso,Group Lasso和嶺回歸都可以看作是本文所提方法的特例。另外，本文的方法也可以用來處理Overlapping Group Lasso問題。

2 數值模擬

例1 (Ω分塊對角)p=100，s*=15，真實的系數向量β*=(0.3,0.3,…,0.3,0,0,…,0)T，預測變量按如下方式生成：

Xj=Z1+0.4εj,Z1～N(0,1), 1≤j≤5Xj=Z2+0.4εj,Z2～N(0,1), 6≤j≤10Xj=Z3+0.4εj,Z3～N(0,1), 11≤j≤15Xj～N(0,1), 16≤j≤100

例2 (Ω帶狀)p=100，β*與例1相同，預測變量(X1,X2,…,Xp)T～N(0,Σ)，其中Σij=0.5|i-j|。此時，若|i-j|=1，有ωii=1.333，ωij=-0.677；若|i-j|>1，則ωij=0。

例3 (Ω稀疏)p=100，預測變量(X1,X2,…,Xp)T～N(0,Ω-1)，其中Ω-1=B+δI。B中非對角線元素以概率0.05取0.5，以概率0.95取0，對角線元素為0。選取δ使得Ω的條件數為p。令β*=Ωγ，其中γ=(γ1,γ2,…,γp)T，且γi=0.1,i=1,2,3,4，否則，γi=0。

表1—表3給出了不同方法在不同的圖結構中的表現，其中l2距離和預測誤差衡量的是模型的估計和預測能力，RFPR和RFNR衡量是模型的變量選擇能力。LG-O和LG表示本文所提方法分別利用真實圖結構和估計圖結構得到的結果。表1—表3中的結果表明：當圖結構比較簡單時，本文方法與其他方法相當(表1所示)；當圖結構比較復雜時，本文方法無論從估計和預測方面還是從模型選擇方面都明顯優于其他方法，與具有Oracle性質的SCAD和MCP方法比較仍具有優越性(表2，表3所示)。

表4給出了不同方法RNMR和RZMR的表現。表4結果表明:即使圖結構比較簡單時，其他方法的RNMR和RZMR也很差。本文的方法在圖結構較簡單時明顯優于其他方法，當圖結構較復雜時，也要優于其他方法，這表明本文方法有效地利用了預測變量圖結構的信息。

表1 例1中不同方法在模型估計、預測和模型選擇能力的比較

表2 例2中不同方法在模型估計、預測和模型選擇能力的比較

表3 例3中不同方法在模型估計、預測和模型選擇能力的比較

表4 不同模型的NMR和ZMR的比較

*注：…表示不可用值，沒有用的預測變量之間沒有相連的邊。

以上模擬結果表明：即使圖結構未知，本文方法仍然能夠有效地利用預測變量圖結構信息，從而提高模型在估計、預測以及變量選擇等方面的表現。

3 實例分析與應用

本節將上述方法應用于公開的乳腺癌實際數據，包括133名受試者22 283個基因表達水平，其中34名為病理完全反應(pCR)受試者，99名為殘留病(RD)受試者。該數據可以通過網址http://Bioinformatics. mdanderson.org/pubdata.html下載。考慮到迭代算法運行的速度以及計算機的限制，將主要比較本文的方法與Lasso方法、嶺回歸方法、Adaptive Lasso方法以及Elastic Net方法的效果。

為了估計精度矩陣，將數據集和測試集分別隨機劃分為大小為112和21的兩部分，然后將整個過程重復100次。每次都采用分層抽樣的方法從對應的組中選取5個pCR個體和16個RD個體作為測試集(大體相當于每組的1/6),其余的個體作為訓練集。在每個訓練集上，利用兩樣本t檢驗選取最顯著的113個基因作為預測變量。注意到，訓練集的樣本量為n=112，比變量的維度p=113稍小，這一點可以用來檢驗p>n時各種方法的表現。利用訓練集估計精度矩陣Ω，基于該精度矩陣利用Graphical Lasso估計預測變量的圖結構G，其圖結構如圖1所示，圖中包含113個節點，190條邊。訓練集用來擬合模型。利用測試集數據計算均方誤差值(FMSE)來評估模型,利用10折交叉驗證選取調節參數[7,17]。

實例分析結果如圖2所示，圖2給出了不同方法在對乳腺癌數據進行建模分析時的平均均方誤差箱線圖。結果可以看出本文的方法在MSE的表現上優于Lasso方法、Adaptive Lasso和Elastic Net方法，比嶺回歸的結果稍差。實際數據分析結果與數值模擬結果一致，說明該方法有效。

圖1 乳腺癌數據的圖結構

圖2 乳腺癌數據集中不同方法的FMSE

4 結論與討論

提出一種基于預測變量圖結構的高維邏輯回歸模型，該模型可以用來對高維圖結構數據或者重疊組結構數據進行邏輯回歸建模，即使預測變量的圖結構未知，本研究的方法仍然能夠利用這種結構提高模型在估計、預測以及變量選擇等方面的表現。另外，目前流行的方法，如Lasso方法、Group Lasso和嶺回歸方法等都可以看作是本文模型的特例。數值模擬和實例分析應用表明：該模型在有限樣本情形下是有效的，并且可廣泛應用于高通量基因數據中存在圖結構的疾病分類研究中。