一種噪聲容錯弱監督矩陣補全的生存分析方法

陳 蕾 邵 楷 林騰濤 陳興國

作為現代統計學的一個重要分支,生存分析旨在建模某個感興趣事件的發生時間.這些感興趣事件通常包括臨床治療中患者的生存時間[1]、機械系統中故障的發生時間[2]以及客戶行為分析中用戶的購買時間[3]等.當前,訓練數據缺乏和數據質量不高是生存分析研究面臨的兩個重要挑戰.一方面,由于數據采集的代價高昂,收集到的樣本數量往往偏少(這一點在醫學領域表現尤為突出),如果再考慮到一些應用領域實例所具有的內在高維特性,那么生存分析問題往往就是一個典型的高維小樣本問題.另一方面,由于觀測周期有限以及研究對象失訪等原因,收集到的數據中還不可避免地存在一些刪失實例,也就是說這些實例的感興趣事件沒有在觀測周期內發生或者由于跟蹤軌跡失效等原因未被觀測到.特別地,為描述方便,對于那些已經觀察到感興趣事件發生的實例,本文將其感興趣事件的發生時間稱為生存時間,相應的實例稱為未刪失實例.而對于其他未觀察到感興趣事件發生的實例,本文將該實例所耗的觀測時間稱為刪失時間,相應的實例稱為刪失實例.通常右刪失是實際應用中最常見的情形,即刪失實例的實際生存時間大于或者等于刪失時間.因此,不失一般性,本文將生存分析的研究對象限定為右刪失實例.

傳統的回歸分析模型通常只能將未刪失實例作為訓練數據,而未刪失實例數量的明顯不足很容易導致模型的過擬合,而另一方面刪失實例的合理利用則有助于改善生存分析模型的泛化性能.因此,一般不采用傳統的回歸模型來處理生存分析問題.Cox 比例風險模型[4]和參數刪失回歸模型[5]是生存分析中兩類應用最為廣泛的模型.但是,這兩類方法均存在較為嚴重的缺陷,所以預測效果并不是很理想.具體地,Cox 模型是建立在比例風險假設上的生存分析模型,其假設實例之間的風險比是常數.因此導致所有實例的生存曲線將呈現相似的形狀.顯然,這種假設在實際應用中是過于嚴格的,并且如果數據中存在具有相同生存時間的實例時,Cox模型還需使用一些近似的方法來處理數據,這有可能會帶來導入偏差(Inducing bias).另一方面,參數刪失回歸模型的預測性能高度依賴于生存時間分布假設的選取.然而,在實際應用中,存在很多影響感興趣事件發生的因素,因此很難選擇出一個合適的理論分布.

近年來,隨著機器學習理論在各個領域的出色表現,研究人員開始引入機器學習方法來研究生存分析問題.相較于傳統的生存分析模型,機器學習能夠從有限的數據中學習到更多的信息,例如數據特征的分布規律及實例特征的抽象表示,同時還具有出色的函數擬合能力.尤其是機器學習中的多任務學習方法[6],能夠幫助模型學習到多個任務之間的共享判別特征,從而提高模型的泛化能力,降低新實例的預測錯誤率.例如,Li 等[7]創新性地將生存分析問題建模成預測各個時間間隔生存狀態的多任務問題,并通過對所設計的線性預測器施加?2,1范數正則化,不僅可以篩選出具有跨任務判別能力的共享特征,還可以緩解高維小樣本問題所帶來的過擬合缺陷,他們所提出的MTLSA/ MTLSA.V2模型在生存分析問題上取得了很好的效果.但是,這些多任務學習方法依然存在一些缺陷,即這些模型均未考慮到數據中可能存在的復雜噪聲所帶來的影響.

為了克服上述缺陷,本文提出一類噪聲容錯弱監督直推式矩陣補全模型(Weakly supervised transductive matrix completion,WSTMC)來預測刪失實例和新實例的生存時間.具體地,基于實例特征潛在的低秩屬性,本文首先將生存分析問題建模為一類多任務直推式矩陣補全(Multi-task transductive matrix completion,MTMC)模型[8].受益于直推式學習機制,該模型不僅可以利用刪失實例作為有限訓練樣本(非刪失實例)的有效補充,而且可以在訓練階段同時探索測試樣本和訓練樣本的特征分布,從而提高模型在測試樣本上的泛化能力.其次,為了克服MTMC 模型的噪聲敏感性,本文引入混合高斯分布(Mixture of Gaussians distribution,MoG)來擬合實際應用中的未知噪聲類型,其動機在于MoG 理論上可以擬合任意連續分布[9].進一步,為了緩解實際應用中高維小樣本所帶來的過擬合缺陷,我們設計了一類新穎的多任務直推式特征選擇機制,以期所提出的模型能在去噪后的特征空間自適應地選擇出跨任務的共享判別特征,從而進一步增強模型的泛化能力.同時,我們還引入了相鄰時間間隔生存狀態的先驗時序穩定性來指導模型生成軟類別標記.此外,我們也設計了一類擬期望最大化迭代優化算法來求解所提出的WSTMC 模型.針對該模型所涉及的多個超參數,我們采用了貝葉斯優化方法來自適應地進行選擇.最后,5 個真實數據集上的實驗結果驗證了所提出的WSTMC 模型優于當前廣泛使用的4 類18 種生存分析方法.

本文結構安排如下:第1 節對相關工作進行討論;第2 節介紹一些數學基礎及矩陣補全預備知識;第3 節為本文主要部分,詳細闡述了所提出的噪聲容錯弱監督直推式矩陣補全(WSTMC)模型及相應的優化算法;第4 節是實驗部分,首先介紹了實驗數據集和所采用的性能評價指標,接著闡述了實驗方法和對比模型,最后分析了實驗結果;第5 節對本文研究內容進行了總結與展望.

1 相關工作

現有的生存分析方法大致可分為如下幾類:1) Cox 比例風險模型;2)參數刪失回歸模型;3)線性模型;4)多任務學習模型.其中,Cox 比例風險模型是生存分析中應用最為廣泛的模型之一.在統計和數據挖掘領域,該模型吸引了眾多研究人員的關注與興趣.近年來,由于各種數據采集和數據分析技術的發展,高維數據在各類實際應用包括生存分析領域頻繁出現,例如醫學領域中微陣列基因表達數據往往包含數千維特征.為了緩解高維數據帶來的過擬合問題,研究人員在傳統的Cox 模型中加入各種正則化項來幫助模型選擇合適的特征以降低特征維度,從而提出了LASSO-COX (Least absolute shrinkage and selection operator Cox)模型[10]、EN-COX (Elastic-net Cox)模型[11]和KEN-COX(Kernel elastic-net Cox)模型[12].它們分別在傳統模型中加入了?1范數正則化項、彈性網正則化項和核化的彈性網正則化項.這種方式擴展了Cox 比例風險模型的應用場景,但并沒有解決該模型中比例風險假設在實際應用中過于嚴苛而可能帶來的欠擬合問題.參數刪失回歸模型是生存分析中另一類常用的重要模型[5,13-15].該模型假設實例的生存時間或者對數生存時間滿足某種特殊的數據分布.其中較常見的假設包括指數分布、威布爾分布(Weibull distribution)、對數(Logistic)分布和極值分布[13]等.此外,一些線性模型還進一步假設生存時間或對數生存時間與特征之間存在線性關聯關系.這類線性模型可以被看成一類特殊的參數刪失回歸模型.例如,結合高斯分布的Tobit[14]模型以及引入Kaplan-Meier 估計的Buckley-James 回歸模型[15].在實際應用中,如何選擇出合適的分布限制了參數刪失回歸模型進一步提高生存分析的預測效果.

近年來,為了松弛傳統統計模型中所依賴的嚴苛假設,即比例風險假設和生存時間分布假設等,一些研究人員開始采用機器學習領域的多任務學習模型來研究生存分析問題.例如,Li 等[7]提出了一種?2,1范數正則化的多任務學習模型MTLSA/MTLSA.V2,該模型利用?2,1范數來共享任務之間的相關性,同時解決高維小樣本數據帶來的過擬合問題;此外還在模型中加入了相鄰時間間隔內生存狀態穩定的先驗條件來約束模型生成標記.但是,該方法僅僅預設數據分布包含單一類型的高斯噪聲,忽略了特征矩陣中可能存在的未知復雜噪聲類型,因此存在一定程度的噪聲敏感性缺陷.受現有模型尤其是近年來多任務學習生存分析模型的啟發,本文提出了一種噪聲容錯的弱監督多任務學習(WSTMC)模型來處理生存分析問題.相較于以往的模型,所提出的WSTMC 模型不僅有效利用了特征矩陣的內在低秩屬性,還能夠通過所引入的MoG 噪聲模型來平滑數據分布中存在的復雜噪聲.此外,不同于傳統的多任務學習方法,我們還在去噪后的特征空間引入了多任務直推式特征選擇機制,幫助模型選擇出跨任務共享的判別特征,這在一定程度上有利于 緩解高維小樣本數據固有的過擬合問題.

2 預備知識

2.1 數學基礎

定義1.瘦型奇異值分解[16].設矩陣X是一個秩為r的m×d維矩陣,則存在兩個單位正交矩陣U ∈Rm×r和V∈Rd×r, 以及對角陣 Σ=diag{σi|1≤i ≤r},其中奇異值σ1≥σ2≥···≥σr＞0 且滿足

式(1)稱為矩陣X的瘦型奇異值分解.

定義2.矩陣范數[16].設秩為r的矩陣X ∈Rm×d存在如式(1)所示的瘦型奇異值分解,則

1)矩陣X的核范數定義為

2)矩陣X的Frobenius 范數定義為

3)矩陣X的?2,1范數定義為

定義3.近鄰算子[17].設F(X)為矩陣X的實值凸函數,則對于任意矩陣M及τ＞0, 函數F(X)的近鄰算子定義為

定義4.閾值收縮算子[17].假設τ＞0, 閾值收縮算子 proxτ‖X‖(M)定義為

其中,Sτ(M)定義為

其中,sign(·)為符號函數.

定義5.核范數近鄰算子[17].設秩為r的矩陣M的奇異值分解為UΣVT, 對任意的τ＞0, 核范數近鄰算子 proxτ‖X‖*(M)定義為

定義6.?2,1范數近鄰算子[18].對于矩陣M ∈Rm×d和任意的τ＞0,其相應的?2,1范數近鄰算子proxτ‖X‖2,1(M)定義為

其中,Jτ(M)定義為

其中,i=1,2,···,m,M(i)表示矩陣M的第i行,‖·‖2表示向量的?2范數.

定理1.近鄰前向后向分裂(Proximal forward backward splitting,PFBS)算法[18].假設F1,F2是兩個下半連續的凸函數,F2在 Rm×d中可微且對某個常數β＞0具有β-Lipschitz 連續梯度,即||?F2(U)-?F2(V)||F≤β||U-V||F, 則凸優化問題

有如下性質:

1) 如果F1+F2是強制的,即lim‖X‖F→+∞(F1(X)+F2(X))=+∞,則該凸優化問題至少存在一個解;

2) 如果F1+F2是嚴格凸的,則該凸優化問題至多存在一個解;

3) 如果F1和F2滿足上述兩個條件,則凸優化問題存在唯一解,且對于任意的初始值X0以及0＜δ ＜2/β用如下方法生成的迭代序列Xk+1收斂到凸優化問題的唯一解

2.2 矩陣補全

矩陣補全可視為壓縮感知理論[19-21]從一維向量空間向二維矩陣空間的一種自然推廣,旨在研究如何在數據不完整的情形下對數據的缺失信息進行填補.標準的矩陣補全問題可以建模為如下形式的秩最小化約束優化模型[22]:

其中,X為待補全的目標矩陣,M為部分元素已知的采樣矩陣,Ω 表示采樣元素的索引集合,PΩ(·)為正交投影算子,表示當 (i,j)∈Ω時,Mij為采樣元素,即

Fazel[23]證明了矩陣核范數是秩函數在矩陣譜范數意義下單位球上的最佳凸逼近.因此,類似于壓縮感知理論中常用的將向量?0范數松弛為向量?1范數的技巧,為了使標準的矩陣補全問題易于求解,一個自然的想法就是利用凸核范數代替非凸秩函數,也就是將原先的秩最小化問題松弛為如下形式的核范數最小化模型[22]:

然而,盡管標準的矩陣補全理論在數據表示、信息重建以及圖像恢復等領域取得了較大成功,但是仍然無法解決弱監督情形下的多標記分類問題.為此,Goldberg 等[8]于2010 年提出了一類新穎的多任務直推式矩陣補全(MTMC)模型并將其成功應用于多標記圖像分類問題.這類MTMC 模型主要建立在兩個假設之上.首先,假設包含m個d維樣本的特征矩陣X和t個任務的任務矩陣Y ∈{0,1,?}m×t之間存在線性依賴關系,也就是說,存在隱含的權重矩陣W0∈Rd×t使得Y=XW0. 其次,假設特征矩陣X滿足內在低秩屬性.基于上述假設,容易推斷出特征-任務堆疊矩陣 [X,Y] 也滿足低秩性.因此標準的MTMC 模型可建模如下:

其中,Z=[ZX,ZY] 表示待求解的特征-任務堆疊矩陣,ZX表示待求解的去噪特征子矩陣,ZY表示待求解的軟標記子矩陣,ΩY表示任務矩陣Y中已知標記的索引集合,Cy(·)表示標記損失函數(通常為logistic 損失函數或平方損失函數).

3 面向生存分析的WSTMC模型

3.1 問題闡述

在生存分析問題中,對于第i個實例,如果該實例未刪失,則可以觀察到它的生存時間S_time(i),反之,可以得到該實例的刪失時間C_time(i).為了方便表達,定義觀察時間O_time(i)如下:

其中,δi∈{0,1}表示刪失狀態指示符,δi=1 表示該實例為刪失實例,δi=0 表示該實例為未刪失實例.通常,可以采用三元組 (xi,O_time(i),δi)來表示生存分析問題中的實例,其中下標i表示實例編號.

本文目標是通過合理建模特征xi和生存時間S_time(i)之間的關系,從而能夠依據實例特征,準確預測該實例的生存時間.然而,由于存在刪失數據,傳統的分類和回歸方法并不適合生存分析問題.受Li 等[7]的啟發,為了更好地利用刪失實例,本文將生存分析問題建模為多任務學習模型.具體地,先將連續的時間進行離散化,也就是將整個研究周期離散化為若干個相等的時間間隔,通過考察各個時間間隔上實例的生存狀態從而間接得到整個研究周期內實例的生存狀態.由此每個實例(xi,O_time(i),δi) 可轉換為三元組 (xi,O_num(i),δi),其中O_num(i)表示第i個實例的觀察時間所對應的時間間隔數.此外,第i個實例在第j個時間間隔上的生存狀態用Yij∈Y來表示.如果在該時間間隔內感興趣事件未發生則Yij=1, 已發生則Yij=0, 狀態未知則使用“?”表示.通過上述方式,生存分析問題轉化為了多任務學習問題,預測每個時間間隔上的生存狀態就是一個學習任務.值得注意的是,對于每個數據集選擇最大的O_num(i)作為多任務學習的任務總數t.圖1 圖示了如何將生存分析問題表述為弱監督多任務學習問題.從圖1 可以看出,對于第i個實例,如果實例是未刪失的,則任務矩陣Y中相應行的生存狀態從開始到第O_num(i)個時間間隔都記為“1”,之后都記為“0”.如果實例是刪失的,則任務矩陣Y中相應行的生存狀態從開始到第O_num(i)個時間間隔都記為“1”,之后都記為“?”.

圖1 生存分析問題建模為弱監督多任務學習問題的圖示Fig.1 Illustration of formulating the survival analysis problem as a weakly supervised multi-task learning problem

綜上所述,可以通過先將生存分析問題闡述為弱監督多任務問題,然后引入交叉熵損失函數將其建模為如下形式的改進MTMC 模型

3.2 建模

然而,直接將生存分析問題闡述為上述多任務直推式矩陣補全(MTMC)問題仍然存在以下幾點不足:1)噪聲容錯性差.MTMC 模型采用Frobenius 范數擬合特征噪聲,已有研究表明Frobenius 范數是高斯分布的最佳逼近[21],然而生存分析數據所面臨的噪聲通常并不是單一的高斯分布且其實際分布類型也往往是未知的[24].2)高維小樣本過擬合.在實際的生存分析問題中,由于研究代價高昂和觀測時間的限制,能收集到的觀測實例往往數量偏少.而受益于數據收集和檢測技術的發展,所能獲取的實例特征則越來越豐富,這加劇了高維小樣本問題所帶來的模型過擬合風險.3)先驗信息融合性不強.MTMC 模型忽略了數據固有的先驗信息,在生存分析問題中,相鄰時間間隔內生存狀態存在時序穩定性,這種先驗時序信息的合理融合通常能有效提升模型的預測性能.

為此,針對噪聲容錯性差的缺陷,我們考慮引入MoG 模型來擬合一般的未知復雜噪聲,其動機來源于已有研究表明MoG 分布能在理論上逼近任意連續分布,這種噪聲建模思想已經在一些典型的機器學習及計算機視覺任務中得到了很好的應用[24-25].具體地,不妨假設特征矩陣中的每一個元素均由兩部分組成,即

其中,Xij,分別表示第i個樣本的第j個觀測特征、真實潛在特征及相應的噪聲.假設噪聲Eij服從獨立的同MoG 分布,即

進一步可得到log 似然函數

我們的目標是通過最大化該log 似然函數解析出特征矩陣X中存在的未知復雜噪聲.為此,結合式(18),易得到如下所示的噪聲容錯MTMC 模型

其次,針對實際應用中高維小樣本所帶來的過擬合缺陷,我們設計了一種新穎的多任務直推式特征選擇機制來自適應地選擇跨任務之間的共享判別特征,從而一定程度上降低原始高維樣本的特征維度,進而緩解模型的過擬合缺陷.為此,首先考慮將MTMC 模型中關于實例特征與標記之間的隱式線性依賴假設ZY=ZXW以正則化項的形式顯式加入目標函數(23),即

其中,W∈Rd×t表示顯式的線性預測器.然后,繼續考慮在上述模型(24)的基礎上引入如下的多任務直推式特征選擇機制

其中,||W||2,1項用來約束線性預測器W保持行稀疏,從而學習到跨任務之間的共享判別特征.特別地,我們注意到這里采用的特征選擇明顯區別于一般的特征選擇機制,一方面我們是在去噪特征空間中實施特征選擇,另一方面所有的實例(包括刪失實例、未刪失實例以及測試實例)均參與了特征選擇.

最后,針對先驗信息融合性不強的缺陷,進一步考慮引入如下的Toeplitz 矩陣S∈Rt×(t-1)并以生成正則化項的形式來誘導軟標記矩陣ZY的

為此,最終構建出噪聲容錯弱監督直推式矩陣補全(WSTMC)模型為

3.3 優化

由于所提出的WSTMC 模型涉及到最大似然估計項L(E,π,σ)的求解,受到期望最大化方法(Experimental maximization,EM)的啟發,我們也考慮引入隱含變量來表征噪聲Eij是否屬于混合高斯分布的第k個分量,從而得到完全數據的log 似然函數

然后,采用EM 方法的思想求解所提出的WSTMC 模型,為此,我們令θ=(Z,W,E,π,σ),首先通過E 步求出隱含變量hijk的期望rijk(這里rijk表示在當前模型參數下第 (i,j)個觀測數據Eij來自第k個分量的概率,也即分量k對觀測數據Eij的響應度):

E-step (求解rijk).

此時，可采用交替求解方法,因此有:

M-step1 (求解πk).

M-step2 (求解).

M-step3 (求解Z,W).

具體地,對于變量Z,根據定理1 可由如下方式求解:

M-step3.1 (求解Z).

對于變量W,同樣根據定理1 可由如下方式求解:

M-step3.2 (求解W).

基于上述分析,所提出的求解WSTMC 模型的擬期望最大化優化算法可概述為算法1.

3.4 時間復雜度分析

求解WSTMC 模型的算法1 是一個迭代優化算法,迭代過程由E-step 和M-step 組成.假設算法1 的外循環迭代次數為NEM, 簡單分析可知Estep 更新rijk所需的時間復雜度為 O(Kmd);基于E-step 的計算結果更新M-step 的π,σ所需的時間復雜度均為 O(Kmd);基于r和σ計算中間變量矩陣B所需的時間復雜度為 O(Kmd);接下來,在更新M-step 的Z,W時又涉及到兩個需要迭代求解的子問題(35a)和(35b),此時不妨假設求解這兩個子問題的內部循環迭代次數為NBPL, 那么對于Mstep 的子問題(35a),求解Z的時間復雜度為O(mdt+m(d+t)min{m,(d+t)}),對于M-step 的另一個子問題(35b),求解W的時間復雜度則為O(mdt+md2).綜上可知,算法1 的時間復雜度為O(NEM(Kmd+NBPL(m(d+t)min{m,(d+t)}+mdt+md2))).為簡單起見,考慮到實際問題中樣本的特征維度d通常大于任務個數t,所以求解模型WSTMC的算法時間復雜度可以簡化為O(NEMNBPL(m2d+md2)).在本文第4 節的實驗部分,我們比較了WSTMC 模型和其他6 種同屬于多任務學習范型的生存分析模型,對于這6 種相關模型的詳細介紹及實驗結果參見第4.4 節及第4.5 節.表1 給出了WSTMC 及其他6 種相關模型的時間復雜度(為簡單起見,這里給出的時間復雜度均假設樣本的特征維度d大于任務個數t).從表1 可以看出,相較于其他6 種多任務學習方法,本文所提出的WSTMC模型具有最高的算法時間復雜度.然而,實驗中發現,每一步并不需要求出子問題的精確解,實際上,只需更新Z與W各一次得到子問題的一個近似解,已足以使算法最終獲得與精確求解子問題時相當的模型性能.

表1 WSTMC 及其他相關模型的時間復雜度比較Table 1 Time complexity comparison of the proposed WSTMC and the other related models

4 實驗

在本節中,首先介紹實驗數據集和評價指標,然后介紹所采用的實驗方法和生存分析比較模型,最后報告實驗結果并展開實驗分析.

4.1 實驗數據集

本文使用了5 個公開的癌癥生存分析數據集,具體包括:NSBCD (Norwegian/Stanford Breast Cancer Data)[29]、DBCD (Dutch Breast Cancer Data)[30]、Lung (Gene Expression Profiles of Lung Adenocarcinoma)[31]、DLBCL (Diffuse Large BCell Lymphoma)[32]和 MCL (Mantle Cell Lymphoma)[33]數據集.這些數據集可以分別從http://user.it.uu.se/～liuya 610/download.html 和http://llmpp.nih.gov/MCL/下載.表2 概述了這些數據集的相關信息,其中“#Instances”列表示數據集所包含的實例數(包括刪失實例和未刪失實例),“#Features”列表示相應數據集中實例的特征數,“ #Censored”列表示數據集中所含的刪失實例數,“#Labels”列表示每個數據集對應的任務個數(即所劃分的時間間隔數,其中NSBCD 和Lung 數據集以“月”作為時間間隔單位;DBCD、MCL 和DLBCL 數據集以“年”作為時間間隔單位).此外,為了表明所采用的數據集是否為高維小樣本問題,我們還在“#Ratios”列記錄了每個數據集中實例數與特征數之間的比例,通常認為當樣本個數比特征維數低一個數量級時即為高維小樣本問題,按照這個標準,容易發現表2 中后四個數據集都屬于高維小樣本問題.

表2 實驗所用數據集概述Table 2 Details of datasets used in this study

4.2 評價指標

由于數據集中存在刪失數據,傳統回歸模型常用的評價指標在生存分析中并不適用.類似于文獻[7],本文也選用了C-index 和Weighted average AUC 兩個指標來評估生存分析性能.其中Cindex 指標側重評估模型在所有任務上的整體回歸性能,而Weighted average AUC 指標則注重評估模型在各個任務上的平均分類性能.兩個指標的具體定義如下:

1) C-index:該指標旨在通過考慮不同事件的相對風險來評估預測值和實際值之間的差異.例如,考慮一對二元組變量()和(),其中表示第i個實例的預測存活時間,Oi表示第i個實例的觀察時間.首先定義一致性概率為

根據定義,對于可以直接預測生存時間的模型來說,C-index 指標計算如下:

其中,I(·)為指示函數.對于多任務類型的生存分析模型,可以通過判斷事件是否發生(以設定閾值的方式)從而計算出生存時間.由于閾值的選取具有主觀性,因此根據Li 等[7]的建議,多任務模型中Cindex 可以通過如下方式計算:

其中,Si表示第i個實例在所有任務上評分之和,即

2) Weighted average AUC:該指標旨在評估模型的整體分類性能,即評估模型能否準確預測出實例在某時間間隔上的生存狀態.將1 當做事件未發生的標記,0 當做事件已發生的標記,則每個時間段就可以看成一個分類任務,整個任務矩陣就相當于t組分類任務.具體地,Weighted average AUC指標可以定義如下:

其中,AUC(i)表示第i個任務的AUC 值,表示第i個任務上已知生存狀態的實例數目.

4.3 實驗方法

為了確定所提出的WSTMC 模型的最佳參數(γ,β,λ,μ,α,c,K),本文采用交叉驗證的方式來評估實驗結果.具體來說,對于實例數目超過150 的生存數據集,我們采用5 折交叉驗證,其他的數據集則采用3 折交叉驗證.為了公平起見,其他對比方法均采用上述的交叉驗證方式.

另一方面,傳統上一般使用網格搜索來尋找合適的參數,但對于本文所提出的WSTMC 模型來說,這將非常耗時.因此,需要一種更加有效的超參數選擇策略.本文采用了基于貝葉斯優化[34]的自動調參方法來選擇最佳的超參數.該方法可以幫助模型選擇出具有更高概率提升預測效果的超參數.具體來說,貝葉斯優化假設存在一個未知的函數ψ=CV_WSTMC(P),其中P∈R7表示一個超參數向量,ψ表示WSTMC 模型對于每組輸入的參數向量進行交叉驗證后取得的預測結果.顯然,最終目的是尋找這個函數的最大值以及對應的超參數向量.首先,貝葉斯優化假設未知函數的結果是通過高斯過程(Gaussian processes,GP)采樣得到的.因此,可以基于歷史記錄(先前的超參數值及其相應的交叉驗證的預測結果)計算出超參數的后驗概率分布.然后,重復迭代下面3 個步驟,直到滿足停止標準(迭代到最大次數,或者連續10 次迭代ψ值都無法得到提升).

步驟1.根據最大化采集函數選擇出下一個最具有“潛力”的評估點(超參數);

步驟2.根據選擇出的參數評估未知函數的預測值;

步驟3.將當前的記錄加入到歷史記錄中,并更新高斯過程.值得注意的是,本文沒有選擇PI(Probability of improvement)作為采集函數,而是使用EI (Expected improvement)作為采集函數來選擇下一次試驗的超參數.雖然PI 策略能夠選擇相對當前預測結果提升概率最大的評估點,但它僅僅反映了提升的概率,并沒有反映提升量的大小.與之相比,EI 策略兩者都能考慮到,還能進一步處理局部和全局之間的關系[35].此外,協方差函數也存在很多種選擇方案[34].本文采用了自動相關性確定(Automatic relevance determination,ARD)作為高斯過程的協方差函數[36],它有助于有效地識別并去除參數向量中不相關的維度.

4.4 對比模型

為了驗證所提出的WSTMC 模型的有效性,我們將其與18 種廣泛使用的生存分析模型進行了比較.這些模型可歸納為如下4 類:基于Cox 的模型、參數刪失回歸模型、線性模型和基于多任務的模型.表3 概述了這18 種比較模型和所提出的WSTMC 模型的相關特點,主要從噪聲容錯機制、直推式學習機制、時序穩定性機制、特征選擇機制和多任務學習機制等5 個方面進行了比較.具體地,對于Cox 模型,我們選擇了傳統Cox 模型,LASSOCOX 模型[10],EN-Cox 模型[11],Cox-?2,1模型和Cox-Trace 模型[37].對于參數模型,選擇了4 種基于不同分布假設的模型,即Weibull、Logistic、Log-logistic 和Log-Gaussian.此外,還與3 種線性模型進行了比較,包括普通的最小二乘(Ordinary least square,OLS)模型[38]、Tobit 模型[14]和RWRSS(Regularized weighted residual sum-of-squares) 模型[35].對于基于多任務的模型,選擇了Multi-LASSO模型、Multi-?2,1模型[27]、MTMC 模型[8]、MTLSA/MTLSA.V2[7]模型和NLMC[28]模型.值得注意的是,OLS 模型、Multi-LASSO 模型和Multi-?2,1模型均無法處理刪失實例.因此,它們只能基于未刪失實例進行學習.為公平起見,比較模型所涉及的所有參數均采用了和WSTMC 模型相同的交叉驗證 策略和參數調整方法.

表3 對比模型的特征比較Table 3 Comparison of characteristics for the competing models

4.5 實驗結果

首先,采用C-index 指標來評估生存分析模型整體回歸性能.表4 報告了不同方法在5 種真實數據集上的實驗結果.對于每種數據集,我們用粗體將最優C-index 值突出顯示,也同時用粗體加下劃線標示出第二優和第三優的C-index 值.從表4 可以看到,WSTMC 模型在所有數據集上均取得了最好的結果.此外,分析實驗結果還可以發現排名前三的結果大多是由基于多任務學習的模型所獲得,這表明通過將原始的生存分析問題轉化為預測各個時間間隔內生存狀態的多任務學習問題,可以有效利用任務之間共享的判別特征來提高預測性能.同時,我們也注意到MTLSA 和MTLSA.V2 模型的性能明顯優于Multi-Lasso 和Multi-l2,1,這反映出能夠在訓練階段利用所有實例(刪失實例和未刪失實例)的模型比只能夠利用單一的未刪失實例的模型效果更好.除此之外,對比WSTMC 模型以及MTMC 和NLMC 模型,由于后兩種矩陣補全模型無法自適應地選擇出跨任務共享的判別特征,也沒能利用時序穩定性這一先驗知識,所以效果始終差于我們所提出的模型.另外,我們也發現,相較于MTLSA 和MTLSA.V2 模型,WSTMC 模型引入混合高斯分布來擬合數據分布中未知的復雜噪聲,有利于消除噪聲對特征選擇的影響,從而進一步降低了預測錯誤率.

接著,我們采用Weighted average AUC 指標來評估生存分析模型在各個任務中的平均分類性能.表5 報告了不同方法在5 種真實數據集上的實驗結果.同樣地,對于每種數據集,最佳結果使用粗體突出顯示,其他前兩名的結果使用粗體加下劃線突出顯示.從表5 可以看出,我們的模型除了在DLBCL 數據集上效果略差外,在其他數據集上始終優于其他對比方法.綜合表4 和表5 的實驗結果,我們可以看到參數刪失回歸模型在兩個評價指標中都沒有取得任何前三名的結果.這表明常用的數據分布無法滿足實際需求,也就是說難以選擇出合適的理論分布來建模實際問題.相比之下,Cox-l2,1模型、MTLSA 模型、MTMC 模型和WSTMC 模型都在兩個評價指標中獲得過前三名.一方面,這驗證了l2,1范數和直推式學習機制有利于緩解高維小樣本帶來的過擬合問題.另一方面,也證實了相較于傳統的生存分析模型,多任務學習方法,尤其是有效地利用了先驗時序穩定機制的多任務學習方法能夠很好地解決生存分析問題.

表4 所提出的WSTMC 模型和其他比較模型在C-index 指標上的性能比較(標準差)Table 4 Comparison of the WSTMC and competing models using C-index (standard deviations)

表5 所提出的WSTMC 模型和其他比較模型在Weighted average AUC 指標上的性能比較(標準差)Table 5 Comparison of the WSTMC and competing models using Weighted average AUC (standard deviations)

4.6 消融實驗

進一步,為了驗證WSTMC 模型中MoG 噪聲容錯機制、多任務直推式特征選擇機制以及時序穩定性機制的有效性,我們也比較了3 種WSTMC的消融模型,即:未引入MoG 噪聲容錯機制的WSTMC-nM 模型、未引入多任務直推式特征選擇機制的WSTMC-nF 模型和未引入時序穩定性機制的WSTMC-nT 模型.具體地,我們通過禁用WSTMC 模型中MoG 噪聲容錯機制、多任務直推式特征選擇機制及時序穩定性機制中的任意一種來驗證這些機制在模型中所起到的實際效果,表6 報告了兩種評價指標C-index 和Weighted average AUC 上5 種相關模型MTMC、WSTMC-nM、WSTMC-nT、WSTMC-nF 以及WSTMC 的性能.實驗結果表明:1)相比于傳統的MTMC 模型,三種消融模型WSTMC-nM、WSTMC-nF 和WSTMC-nT 在兩種評價指標C-index 和Weighted average AUC 上都無一例外地優于MTMC,這表明引入MoG 噪聲容錯機制、多任務直推式特征選擇機制以及時序穩定性機制是有助于提高傳統MTMC 模型性能的;2)相比于同時引入3 種機制的WSTMC 模型,消融模型WSTMC-nM、WSTMCnF 和WSTMC-nT 在兩種評價指標C-index 和Weighted average AUC 上性能均弱于WSTMC,進一步驗證了3 種機制的同時引入優于任意兩種機制的組合引入,從而實驗上表明我們所提出的WSTMC 模型中3 種機制的引入是合理有效的.

表6 在兩種評價指標C-index 和Weighted average AUC 上的消融性實驗性能比較(標準差)Table 6 Comparison of the ablation experiments using C-index and Weighted average AUC (standard deviations)

表6 在兩種評價指標C-index 和Weighted average AUC 上的消融性實驗性能比較(標準差) (續表)Table 6 Comparison of the ablation experiments using C-index and Weighted average AUC(standard deviations) (continued)

5 總結與討論

針對生存分析領域通常面臨的高維小樣本和噪聲敏感問題,本文提出了一類新穎的噪聲容錯弱監督直推式矩陣補全(WSTMC)模型.首先將原始的生存分析問題建模成一類傳統的多任務直推式矩陣補全(MTMC)模型,然后引入高斯混合分布來擬合數據中未知的復雜噪聲.同時為了緩解高維小樣本所帶來的過擬合缺陷,我們還設計了一類適用于去噪特征空間的多任務直推式特征選擇機制,以期篩選出跨任務共享的判別特征,從而提高模型的泛化性能.相較于傳統的MTMC 模型,我們提出的WSTMC 模型一定程度上克服了MTMC 模型噪聲容錯性差、泛化性能不足以及先驗信息融合性不強的缺陷.最后,多個真實數據集上的實驗結果證實了所提出的WSTMC 模型優于其他廣泛使用的生存分析方法.

然而,我們的方法也仍然存在改進的空間,比如無論是傳統的MTMC 模型還是我們的WSTMC 模型都是基于實例特征與實例標記的線性依賴假設,這種假設與實際情形未必相符,如果能將我們的模型推廣至非線性假設情形,將極有可能提升生存分析的預測性能.我們將在后續研究中嘗試采用核化的方法來放寬已有模型的線性依賴假設,并同時嘗試引入其他類型的非線性依賴實現機制.