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

陳 蕾 邵 楷 林騰濤 陳興國

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

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

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

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

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

1 相關(guān)工作

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

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

2 預備知識

2.1 數(shù)學基礎(chǔ)

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

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

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

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

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

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

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

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

其中,Sτ(M)定義為

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

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

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

其中,Jτ(M)定義為

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

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

有如下性質(zhì):

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

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

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

2.2 矩陣補全

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

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

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

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

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

3 面向生存分析的WSTMC模型

3.1 問題闡述

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

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

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

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

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

3.2 建模

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

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

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

進一步可得到log 似然函數(shù)

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

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

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

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

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

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

3.3 優(yōu)化

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

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

E-step (求解rijk).

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

M-step1 (求解πk).

M-step2 (求解).

M-step3 (求解Z,W).

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

M-step3.1 (求解Z).

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

M-step3.2 (求解W).

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

3.4 時間復雜度分析

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

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

4 實驗

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

4.1 實驗數(shù)據(jù)集

本文使用了5 個公開的癌癥生存分析數(shù)據(jù)集,具體包括: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]數(shù)據(jù)集.這些數(shù)據(jù)集可以分別從http://user.it.uu.se/～liuya 610/download.html 和http://llmpp.nih.gov/MCL/下載.表2 概述了這些數(shù)據(jù)集的相關(guān)信息,其中“#Instances”列表示數(shù)據(jù)集所包含的實例數(shù)(包括刪失實例和未刪失實例),“#Features”列表示相應數(shù)據(jù)集中實例的特征數(shù),“ #Censored”列表示數(shù)據(jù)集中所含的刪失實例數(shù),“#Labels”列表示每個數(shù)據(jù)集對應的任務(wù)個數(shù)(即所劃分的時間間隔數(shù),其中NSBCD 和Lung 數(shù)據(jù)集以“月”作為時間間隔單位;DBCD、MCL 和DLBCL 數(shù)據(jù)集以“年”作為時間間隔單位).此外,為了表明所采用的數(shù)據(jù)集是否為高維小樣本問題,我們還在“#Ratios”列記錄了每個數(shù)據(jù)集中實例數(shù)與特征數(shù)之間的比例,通常認為當樣本個數(shù)比特征維數(shù)低一個數(shù)量級時即為高維小樣本問題,按照這個標準,容易發(fā)現(xiàn)表2 中后四個數(shù)據(jù)集都屬于高維小樣本問題.

表2 實驗所用數(shù)據(jù)集概述Table 2 Details of datasets used in this study

4.2 評價指標

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

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

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

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

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

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

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

4.3 實驗方法

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

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

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

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

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

4.4 對比模型

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

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

4.5 實驗結(jié)果

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

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

表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 噪聲容錯機制、多任務(wù)直推式特征選擇機制以及時序穩(wěn)定性機制的有效性,我們也比較了3 種WSTMC的消融模型,即:未引入MoG 噪聲容錯機制的WSTMC-nM 模型、未引入多任務(wù)直推式特征選擇機制的WSTMC-nF 模型和未引入時序穩(wěn)定性機制的WSTMC-nT 模型.具體地,我們通過禁用WSTMC 模型中MoG 噪聲容錯機制、多任務(wù)直推式特征選擇機制及時序穩(wěn)定性機制中的任意一種來驗證這些機制在模型中所起到的實際效果,表6 報告了兩種評價指標C-index 和Weighted average AUC 上5 種相關(guān)模型MTMC、WSTMC-nM、WSTMC-nT、WSTMC-nF 以及WSTMC 的性能.實驗結(jié)果表明:1)相比于傳統(tǒng)的MTMC 模型,三種消融模型WSTMC-nM、WSTMC-nF 和WSTMC-nT 在兩種評價指標C-index 和Weighted average AUC 上都無一例外地優(yōu)于MTMC,這表明引入MoG 噪聲容錯機制、多任務(wù)直推式特征選擇機制以及時序穩(wěn)定性機制是有助于提高傳統(tǒng)MTMC 模型性能的;2)相比于同時引入3 種機制的WSTMC 模型,消融模型WSTMC-nM、WSTMCnF 和WSTMC-nT 在兩種評價指標C-index 和Weighted average AUC 上性能均弱于WSTMC,進一步驗證了3 種機制的同時引入優(yōu)于任意兩種機制的組合引入,從而實驗上表明我們所提出的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 上的消融性實驗性能比較(標準差) (續(xù)表)Table 6 Comparison of the ablation experiments using C-index and Weighted average AUC(standard deviations) (continued)

5 總結(jié)與討論

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

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