覃 婷 王 彤

當(dāng)用基因表達數(shù)據(jù)預(yù)測生存情況時，基因數(shù)遠遠超過了樣本例數(shù)。除了高維度以外，基因表達之間通常存在著某種未知的相關(guān)，其增加了解釋變量之間的共線性。基因表達數(shù)據(jù)存在的小樣本、高維度、強相關(guān)的特點給生存預(yù)測帶來了困難。因此根據(jù)基因數(shù)據(jù)做生存預(yù)測時，首先需要對基因表達數(shù)據(jù)進行降維或者調(diào)整，從而更加有效而準確的進行參數(shù)估計。有監(jiān)督的主成分分析(supervised principal component analysis，SuperPC)和偏 Cox回歸(partial least squares Cox regression)是其中的兩種降維方法。本文將通過模擬研究和對國際上公開的三個基因數(shù)據(jù)集進行分析，以探討這兩種方法用于高維數(shù)據(jù)生存預(yù)測表現(xiàn)的優(yōu)劣，為得到更精確的預(yù)后估計和改進治療策略提供依據(jù)。

基本原理與方法

1.基因微陣列數(shù)據(jù)的標識與比例風(fēng)險模型

假設(shè)有一組包含著截尾數(shù)據(jù)的基因微陣列生存數(shù)據(jù)有 n 個個體，(yi，δi，xi);i=1，…，n。其中 yi為個體i的失效時間，yi可以是完全數(shù)據(jù)，也可是截尾數(shù)據(jù);δi是一個指示變量，當(dāng)δi=1時為完全數(shù)據(jù)，而當(dāng)δi=0為截尾數(shù)據(jù);xi=(xi1，…，xip)T為個體i的自變量向量。

令Y為生存時間。生存函數(shù)定義為S(y)=P(Y＞y)，是某個體在時間y時刻依然存活的概率。風(fēng)險函數(shù)測量在y時刻存活的個體，在下個很小的時間段內(nèi)死亡的瞬時風(fēng)險。比例風(fēng)險模型表示為

h(y，X)=h0(y)exp(XTβ) (1)其中h0(y)是一個非指定的基準風(fēng)險函數(shù)。

模型的參數(shù)向量^β通常可以取最大偏對數(shù)似然得到，基準生存函數(shù)H0(y)可用Breslow估計，表示為 ^H0(y)。對于一個表達譜為~X的新樣本，根據(jù)已知的參數(shù)和基線風(fēng)險估計求出其風(fēng)險函數(shù)和生存函數(shù)，

2.有監(jiān)督的主成分分析

SuperPC 是由 Bair和 Tibshirani等人提出的〔1－2〕，它改進常規(guī)主成分分析無法保證所選擇的主成分與病人的生存相關(guān)的缺點，在降維的時候考慮了生存時間，其核心思想就是只對與生存時間密切相關(guān)的基因進行主成分分析。

該方法首先將每一個基因分別代入單變量Cox模型h(t|x)=h0(t)exp(βx)，以檢驗它們對生存的影響。然后對其進行基于偏似然函數(shù)的參數(shù)估計與假設(shè)檢驗，檢驗方法為似然比檢驗。將基因按照檢驗所得到的P值從小到大排序，然后根據(jù)交叉驗證法挑選出前λ1百分比的基因組成一個簡化矩陣Xθ。采用奇異值分解法(singular value decomposition，SVD)對這個簡化矩陣進行主成分分析。

假設(shè)X矩陣的列已經(jīng)被中心化，均數(shù)為0。那么n×p矩陣X的奇異值分解寫作:

其中U是一個n×n的正交陣，V是p×p正交陣，D是一個以奇異值dj為對角元素的n×p對角陣，r=min(n，p)是X的秩，非零奇異值的數(shù)目與X矩陣的秩相等，d1≥d2≥…≥dr＞0。

那么，簡化矩陣Xθ的奇異值分解寫作:

令 Uθ=(uθ，1，uθ，2，…，uθ，r)，稱 uθ，1為 X 的第一有監(jiān)督的主成分，依此類推。如果僅取一個成分，即擬合一個應(yīng)變量為y和自變量為uθ，1的Cox比例風(fēng)險模型，得

從公式(3)變換得到(注意到正交陣V'V=I)，

因為 uθ，1是 Xθ的一個線性組合:uθ，1=Xθwθ，1，所以模型(4)可以看作是一個利用了Xθ中的所有自變量的受限模型:

假如有一個新的基因數(shù)據(jù)集x*，對其進行生存預(yù)測，步驟如下:

3.偏Cox回歸

在基因數(shù)據(jù)的生存預(yù)測方法中，基于偏最小二乘的生存分析是一個重要的家系。利用偏最小二乘方法進行降維，即從原始變量中提取偏最小二乘成分，然后將提取的這些線性成分應(yīng)用于標準的Cox回歸進行生存預(yù)測的方法，就稱之為偏Cox回歸。

偏Cox回歸算法有很多種，這里采用的是Nyg?rd提出的算法〔3〕。該算法主要是通過將生存問題轉(zhuǎn)換為廣義線性回歸問題，然后依照廣義線性模型的迭代再加權(quán)偏最小二乘算法提取PLS成分，從而實現(xiàn)高維數(shù)據(jù)的降維，然后將所得到的參數(shù)估計以及提取的PLS成分代入Cox比例風(fēng)險模型中，進行生存預(yù)測。由于這種算法將基準風(fēng)險增量的估計與PLS降維分開，使得PLS的成分僅為基因表達譜的線性組合，更符合生物解釋。

由于該算法只是對^η進行了部分更新，可解決收斂速度過慢的問題。且分開估計協(xié)方差效應(yīng)和基線風(fēng)險增量避免了數(shù)據(jù)維度的擴張，加上提取的PLS成分數(shù)量小，所以計算速度通常很快，節(jié)約了計算時間。

4.根據(jù)交叉驗證法選擇模型調(diào)整參數(shù)

預(yù)測方法的模型復(fù)雜程度是由估計調(diào)整參數(shù)來決定的。調(diào)整參數(shù)的估計方法有很多種，最經(jīng)常使用的就是交叉驗證。在本文中，我們采用Verweij和van－Houwelingen提出的交叉驗證準則〔4〕，這種準則是建立在Cox偏對數(shù)似然的基礎(chǔ)上的。

首先將數(shù)據(jù)分成等大小的K個部分(1＜K≤n)，每個第i(1≤i≤K)次的交叉驗證都會將第i個層剔除，只用剩下的層來訓(xùn)練模型，根據(jù)訓(xùn)練好的模型來估計被剔除的第i個層的預(yù)測信息，重復(fù)K次，這樣每一個部分都做了并且只做了一次驗證組。令l(β)表示全部數(shù)據(jù)的Cox對數(shù)似然，l(－i)(β)表示剔除第i個層數(shù)據(jù)的對數(shù)似然，將第i個部分對似然的貢獻定義為li(β)=l(β) － l(－i)(β)，使得 l(－i)(β)最大化的 β 估計值表示為β^(－i)。假設(shè)似然成分是獨立的，那么li(β)就簡單的等于第i個部分的貢獻，并且l(β)，K折交叉驗證的最大對數(shù)似然為 CVL =，隨著λ的變化，每個最大對數(shù)似然CVL也在變化，其中最大的CVL所對應(yīng)的λ為最優(yōu)調(diào)整參數(shù)。

在實踐中，通常會給出調(diào)整參數(shù)的范圍，在給定的范圍內(nèi)找最優(yōu)調(diào)整參數(shù)。SuperPC的調(diào)整參數(shù)λ=(λ1，λ2)是雙變量的，λ1代表單變量分析中與生存時間有關(guān)的基因子集占基因全集的百分比，λ2代表將要選擇的主成分數(shù)的范圍。

5.模型評價

一般來說，預(yù)報因子的真正性能評價應(yīng)該用一個獨立的數(shù)據(jù)來完成。但是在缺乏獨立數(shù)據(jù)的情況下，可以通過交叉驗證來進行模型評價，即將每一個數(shù)據(jù)集按2:1的比例隨機分成訓(xùn)練數(shù)據(jù)集和驗證集。訓(xùn)練組樣本的基因表達和生存數(shù)據(jù)被用來構(gòu)建預(yù)測模型，驗證集用來評估模型的性能。為了避免依賴于訓(xùn)練集和驗證集的選擇，需進行重復(fù)分組，以所有評價標準的結(jié)果的中位數(shù)和四分位數(shù)間距來估計預(yù)測模型的性能。本論文主要的評價標準為決定系數(shù)R2。

決定系數(shù)是驗證組中的生存數(shù)據(jù)可以被預(yù)報因子解釋的那部分變異所占的百分比。預(yù)測性能良好的預(yù)報因子可以解釋驗證組生存數(shù)據(jù)的絕大部分變異。在傳統(tǒng)的回歸背景下，R2=1－殘差平方和/總平方和，因此它的取值范圍在0到1之間。然而這個定義在數(shù)據(jù)存在刪失的情況下不能使用，因此，Nagelkerke給出了一個可以用在Cox比例風(fēng)險模型中R2統(tǒng)計量〔5〕

其中l(wèi)(．)表示對數(shù)似然函數(shù)。R2越大，預(yù)報器的預(yù)測性能越好。

模擬分析

為了進一步驗證上述模型的預(yù)測性能，我們根據(jù)基因數(shù)據(jù)的特點設(shè)計了模擬實驗，用上述方法對模擬實驗產(chǎn)生的數(shù)據(jù)進行分析，根據(jù)模型預(yù)測評價標準對它們的預(yù)測性能進行比較。模擬條件設(shè)定如下:

(1)協(xié)變量矩陣X:

生成100×1000的基因協(xié)變量矩陣，每一行表示一個病人，每一列表示一個基因。這些協(xié)變量服從多元正態(tài)分布，均數(shù)向量為0。將數(shù)據(jù)分成十塊等大小的基因塊∑b，令它們的方差協(xié)方差矩陣的對角元素為，非對角元素為。因此∑對應(yīng)于基因表達的類別，這樣不同類別的基因表達是獨立的，但是在同一個類別中的基因表達有同樣的兩兩相關(guān)。在模擬中，我們令ρ分別等于0.3，0.6，0.9以觀察不同相關(guān)程度對結(jié)果的影響，同時評判三種不同方案的方差:(a)=…==1，即所有的基因表達的方差相等，(b)==2，=…==1，即前兩塊的基因表達的變異更大，(c)==1/2=…==1，即前兩塊的基因變異比其他塊的要少。

我們首先產(chǎn)生服從(0，1)均勻分布的隨機數(shù)S，令生存函數(shù)S(t)=S，利用用產(chǎn)生相應(yīng)的生存時間t。

(4)生成刪失指示變量:

產(chǎn)生一組隨機數(shù)，服從二項分布，發(fā)生1的概率為0.8，也就是截尾為20%，為了觀察截尾比例是否對降維方法產(chǎn)生影響，我們改變了截尾比例，截尾比例分別為20%，50%。

按照上述的實驗設(shè)計產(chǎn)生訓(xùn)練組數(shù)據(jù)，然后以同樣的條件產(chǎn)生驗證組數(shù)據(jù)。分別應(yīng)用SuperPC和偏

(2)參數(shù)的設(shè)定:

每個回歸系數(shù)對應(yīng)于它對應(yīng)變量的影響。在本文中，參數(shù)的設(shè)定如下:當(dāng) 1≤j≤100，βj=0.01，當(dāng) 101≤j≤200，βj從0.01 到0.1，每0.01 為一個步長。201≤j≤1000，βj=0，表示在基因矩陣中只有少數(shù)一些相關(guān)的協(xié)變量，大部分的協(xié)變量都是無關(guān)的。

(3)生成生存時間:Cox回歸方法對訓(xùn)練組進行建模，然后用驗證組數(shù)據(jù)進行預(yù)測評價，在最優(yōu)調(diào)整參數(shù)的條件下評價模型的預(yù)測性能，評價標準為R2，從評價標準的中位數(shù)以及離散程度來對模型進行評價。以上過程重復(fù)200次。模擬結(jié)果以箱式圖矩陣表示。

箱式圖矩陣的行代表影響生存的前兩塊基因的方差，列代表3個不同的相關(guān)系數(shù)。spcr為SuperPC方法，pls為偏Cox回歸方法。censored=0.2和censored=0.5分別表示刪失比例為20%和50%的模擬數(shù)據(jù)的結(jié)果。

由圖1可見，根據(jù)模型評價標準，SuperPC的預(yù)測性要優(yōu)于偏Cox;當(dāng)相關(guān)系數(shù)相同時，兩種方法得到的R2都是隨著方差的增大而增大，提示基因的方差越大的時候，基因塊越容易被識別。方差相同時，相關(guān)系數(shù)越大，得到的R2越大。隨著截尾比例的增大，模型的預(yù)測性能會變差，說明模型的預(yù)測能力會受到刪失比例的影響。

實例分析

用上述兩種方法對國際上三個公開的基因數(shù)據(jù)集進行了分析，它們分別是Van't Veer等的乳腺癌數(shù)據(jù)〔6〕，78個病人，4 751個基因，觀察事件的結(jié)局是乳腺癌是否轉(zhuǎn)移，截尾比例為56.4%;Beer等人的肺癌數(shù)據(jù)集〔7〕有86個病人和7 129個基因，觀察事件結(jié)局為死亡，截尾比例為72.6%;Bullinger等人急性髓系白血病(acute myeloid leukemia)的數(shù)據(jù)〔8〕，116 個病人，6 283個基因，觀察事件結(jié)局為死亡，截尾比例為42.24%。

因為對于每個數(shù)據(jù)集，采用不同的評價標準，最佳的預(yù)測方法有可能是不同的〔9〕。因此，對于一個實例數(shù)據(jù)，首先要看用這兩種方法進行分析的預(yù)測性能如何，然后根據(jù)預(yù)測結(jié)果挑選合適的預(yù)測模型。

首先按2:1的比例將數(shù)據(jù)隨機分成訓(xùn)練組和驗證組:訓(xùn)練組用于構(gòu)造模型，而驗證組用來對模型的預(yù)測性能進行評價。為了保證預(yù)測結(jié)果評價的客觀性，避免數(shù)據(jù)任意分割導(dǎo)致的預(yù)測偏差，按上述方法重復(fù)將數(shù)據(jù)集隨機分割200次。結(jié)果見圖2。

圖中bc代表乳腺癌數(shù)據(jù)，lc代表肺癌數(shù)據(jù)，aml代表急性髓系白血病數(shù)據(jù)，spcr為是SuperPC方法的分析結(jié)果，pls為偏Cox方法的分析結(jié)果。

R2值越大，模型可以解釋數(shù)據(jù)的變異部分越大，模型的預(yù)測性能越好。從圖2中我們可以看到，對于乳腺癌數(shù)據(jù)，SuperPC方法的表現(xiàn)要優(yōu)于偏Cox回歸。而對于肺癌數(shù)據(jù)和急性髓系白血病數(shù)據(jù)，則偏Cox回歸的表現(xiàn)要優(yōu)于SuperPC方法。

圖1 模擬數(shù)據(jù)分析結(jié)果矩陣圖

圖2 實際數(shù)據(jù)分析結(jié)果圖

因為我們對這3個數(shù)據(jù)進行了隨機分割(分割成訓(xùn)練集和驗證集)，分割了200次，每次分割得到的最優(yōu)調(diào)整參數(shù)的取值是不一樣的。因此，這里我們給出了三個實例數(shù)據(jù)分析中，兩種方法各自所選擇的最優(yōu)調(diào)整參數(shù)的分位數(shù)的表。

從表1中我們可以看到，SuperPC提取的成分數(shù)大于等于偏Cox回歸，并且提取的成分數(shù)不穩(wěn)定。SuperPC提取的基因子集占原基因集的比例很小，即簡化矩陣遠小于原始矩陣。所以雖然SuperPC引入的成分數(shù)多，但每個成分中包含的自變量數(shù)目要少于全基因集。偏Cox方法提取的成分數(shù)穩(wěn)定，用很少的成分就可以解釋原始變量的大部分變異。

表1 實際數(shù)據(jù)分析中選擇的最優(yōu)調(diào)整參數(shù)的分位數(shù)

根據(jù)上面的模型預(yù)測性能評價，挑選各個數(shù)據(jù)對應(yīng)的最適方法來對整個數(shù)據(jù)集進行了模型擬合就得到了預(yù)測方程?？梢杂嬎泐A(yù)后指數(shù)來估計病人的預(yù)后情況，預(yù)后指數(shù)的公式為PI=X^β。當(dāng)?shù)玫叫碌牟∪说幕驍?shù)據(jù)~X的時候，就可以根據(jù)預(yù)后指數(shù)PI=~X^β來對其分類:按照預(yù)后指數(shù)的中位數(shù)，將病人分成兩組，預(yù)后指數(shù)超過中位數(shù)的，為高危組，可考慮相對積極的治療策略;低于中位數(shù)的，為低危組，可以采用相對保守的治療方案。

討 論

有監(jiān)督的主成分分析和偏Cox回歸都是將Cox比例風(fēng)險模型與降維技術(shù)結(jié)合起來，適用于基因數(shù)據(jù)的生存預(yù)測。這兩種方法都對主成分回歸進行了改進，利用了應(yīng)變量的信息。它們的區(qū)別在于有監(jiān)督的主成分分析是對基因子集進行特征提取，而偏Cox回歸則是對基因的全集進行特征提取。

從模擬實驗中可見:(1)隨著影響生存的基因塊的方差的增大，兩種方法的預(yù)測性能變好，這是因為當(dāng)自變量的變異更大的時候，被提取的信息也更多，可提高回歸估計的精度和穩(wěn)定性，故而當(dāng)影響生存的基因塊的方差大的時候較容易被識別出來。(2)隨著組內(nèi)相關(guān)系數(shù)ρ的增高，兩種方法的預(yù)測性能都變好。因為只有存在著相關(guān)，才能夠把維度降下來，特征值會隨著相關(guān)的增高而增大，那么前幾個主成分的方差就會很大，這樣只需要幾個主成分，就能夠很好的解釋原始變量的變異。隨組內(nèi)相關(guān)系數(shù)的增高，預(yù)測方法性能變好的同時還應(yīng)注意到，當(dāng)基因的相關(guān)度較高的時候，預(yù)測方法對選擇正確的基因比較不敏感，對生存沒有影響卻與對生存有影響的基因高度相關(guān)的基因也可能被納入。(3)隨著刪失比例的增加，兩種方法的預(yù)測性能變差，說明預(yù)測方法會受到刪失比例的影響。

在實例分析中，根據(jù)判斷標準，不同的數(shù)據(jù)集最優(yōu)預(yù)測方法不同。模擬研究和實例分析中，有監(jiān)督的主成分分析提取的成分數(shù)要大于等于偏Cox回歸。然而我們也應(yīng)該注意到，因為這種主成分是有監(jiān)督的，所以雖然引入的成分數(shù)多，但每個成分中包含的自變量數(shù)目要少于全基因集。

在我們的模擬研究中，這兩種方法的預(yù)測性能很好，用我們的模型能夠解釋原始變量的變異比例較大，決定系數(shù)甚至達到了0.9的情況。但是，實例分析中可以看到，決定系數(shù)沒有那么大。造成這種差異的原因可能是因為模擬實驗設(shè)計的時候，基因塊之間設(shè)定了相關(guān)系數(shù)，各個基因塊中的基因兩兩相關(guān)，而基因塊之間不相關(guān)，并且對回歸參數(shù)也進行了設(shè)定，只是前兩塊基因與生存有關(guān)，數(shù)據(jù)生成具有一定的規(guī)律。而實際的微陣列數(shù)據(jù)遠比模擬數(shù)據(jù)要復(fù)雜得多，基因之間以未知的方式相關(guān)，并且微陣列數(shù)據(jù)中的協(xié)變量數(shù)千至數(shù)萬個，簡單的模擬不能夠捕獲這種復(fù)雜的關(guān)系。并且模擬中為了計算的方便，只是模擬了1 000個協(xié)變量，樣本含量固定取100個，這與實際微陣列數(shù)據(jù)的樣本含量和自變量個數(shù)相比，樣本含量的比例要大得多。實例數(shù)據(jù)分析，因為是隨機拆分數(shù)據(jù)，最后得到的訓(xùn)練組數(shù)據(jù)和驗證組數(shù)據(jù)的截尾比例可能會比原來的數(shù)據(jù)要高，這也是影響結(jié)果的一個原因之一。

在模擬研究中，有監(jiān)督的主成分分析的預(yù)測性能要優(yōu)于偏Cox回歸。但本文介紹的偏Cox回歸算法因其在收斂性上進行了改進，從而計算速度很快，比有監(jiān)督的主成分分析方法節(jié)約了很多時間。對于高維數(shù)據(jù)的分析來說，計算時間短也是個很重要的優(yōu)點。

1.Bair E，Tibshirani R.Semi－supervised methods to predict patient survival from gene expression data.PLoS Biology，2004，2:511－522.

2.Bair E，Hastie T，Paul D，et al.Prediction by supervised principal components.Journal of the American Statistical Association，2006，101:119－137.

3.Nyg?rd S，Borgan O，Lingiaerde OC，et al.Partial least squares Cox regression for genome－wide data.Lifetime Data Anal，2008，14:179－195.

4.Verweij PJMvan，Houwelingen HC，Cross－validation in survival analysis.Stat Med，1993，12:2305－2314.

5.Nagelkerke NJS.A note on a general definition of the coefficient of determination.Biometrika，1991，78:691－692.

6.Van't Veer LJ，Dai H，Van de Vijver.Gene expression profiling predicts clinical outcome of breast cancer.Nature，2002，415:530－536.

7.Beer DG，Kardia SL，Huang CC，et al.Gene－expression profiles predict survival of patients with lung adenocarcinoma，2002，Nat Med 8:816－824.8.Lars Bullinger MD，Konstanze D?hner MD，Eric Bair，et al.Use of geneexpression profiling to identify prognostic subclasses in adult acute myeloid leukemia.Massachusetts Medical Society，2004，350 16:1605－1616.

9.W.van Wieringen，D.Kun，R.Hampel，et al.Survival prediction using gene expression data:a review and comparison.Computational Statistics and Data Analysis，2009，53:1590－1603.