中國(guó)醫(yī)學(xué)科學(xué)院基礎(chǔ)醫(yī)學(xué)研究所/北京協(xié)和醫(yī)學(xué)院基礎(chǔ)學(xué)院統(tǒng)計(jì)學(xué)教研室(100005) 王子興 申郁冰 姜晶梅

醫(yī)學(xué)中常常需要評(píng)價(jià)某種標(biāo)志物、影像指標(biāo)或評(píng)分工具對(duì)疾病發(fā)生、診斷和預(yù)后的預(yù)測(cè)性能[1]。受試者工作特征(receiver operating characteristics，ROC)曲線分析是面向上述研究目標(biāo)的一類統(tǒng)計(jì)學(xué)方法，其基于“金標(biāo)準(zhǔn)”將疾病狀態(tài)進(jìn)行二分類，進(jìn)而綜合分析待評(píng)價(jià)指標(biāo)各閾值所對(duì)應(yīng)的靈敏度、特異度信息。然而疾病狀態(tài)與時(shí)間有關(guān)，可經(jīng)歷“從無(wú)到有”的變化過(guò)程[2]。這種依時(shí)間而改變的結(jié)局變量給預(yù)測(cè)性能的評(píng)價(jià)帶來(lái)兩方面問(wèn)題：①基于不同時(shí)間截面的分析會(huì)有不同的疾病狀態(tài)分布；②疾病狀態(tài)可發(fā)生刪失，且刪失比例隨研究時(shí)間的延長(zhǎng)而增高[3-4]。為實(shí)現(xiàn)對(duì)該類含刪失的時(shí)間-狀態(tài)數(shù)據(jù)的預(yù)測(cè)性能評(píng)價(jià)，經(jīng)典ROC曲線和生存分析方法進(jìn)行結(jié)合產(chǎn)生了依時(shí)ROC(time-dependent ROC)分析方法[5]。依時(shí)ROC在過(guò)去20年內(nèi)得到了極為迅速的豐富和發(fā)展，并應(yīng)用于疾病預(yù)防、藥物治療等領(lǐng)域。

本文從非參、半?yún)⒔嵌瓤偨Y(jié)依時(shí)ROC分析的方法學(xué)研究及其對(duì)應(yīng)的軟件實(shí)現(xiàn)方式，以便增進(jìn)讀者對(duì)依時(shí)ROC方法的了解并在臨床轉(zhuǎn)化中的應(yīng)用提供參考。

基本定義

1.病例和對(duì)照的三種定義

依時(shí)ROC除標(biāo)志物X閾值c外，還需結(jié)合時(shí)間閾值t進(jìn)行分析(圖1)，由此產(chǎn)生三種依時(shí)ROC的“病例”和“對(duì)照”定義，這些定義最早由Heagerty和Zheng于2005年總結(jié)[6]并得以沿用。

圖1 依時(shí)ROC三種定義示意圖

累積/動(dòng)態(tài)(cumulative/dynamic，CD)定義下，累積病例為研究起始點(diǎn)到t時(shí)刻期間經(jīng)歷事件的個(gè)體(圖1中的A、B和E)，動(dòng)態(tài)對(duì)照為t時(shí)刻仍未發(fā)生事件的個(gè)體(圖1中C和F)?？梢?jiàn)CD定義下病例和對(duì)照集均隨著時(shí)間變化而改變，某個(gè)體可能在前后分別作為對(duì)照和病例，故存在信息被重復(fù)利用的問(wèn)題[7]。然而CD定義具備鮮明的臨床意義，可用于預(yù)測(cè)個(gè)體是否在指定時(shí)間內(nèi)發(fā)生感興趣結(jié)局，故得到了最大程度的方法學(xué)發(fā)展和應(yīng)用[2]。

時(shí)點(diǎn)/動(dòng)態(tài)(incident/dynamic，ID)定義下，時(shí)點(diǎn)病例為指定時(shí)點(diǎn)t上恰好發(fā)生事件的個(gè)體(圖1中的A)，而動(dòng)態(tài)對(duì)照同前。ID定義可利用其與生存分析中風(fēng)險(xiǎn)集概念的對(duì)應(yīng)關(guān)系得到一些簡(jiǎn)化估計(jì)方法。

時(shí)點(diǎn)/靜態(tài)(incident/static，IS)定義下，時(shí)點(diǎn)病例同前，而靜態(tài)對(duì)照為在某固定隨訪期(O，t*)內(nèi)未發(fā)生事件的個(gè)體(圖1中的C)，其中t*需足夠長(zhǎng)來(lái)觀察終點(diǎn)事件。IS定義主要面向篩查研究。

2.符號(hào)定義

以n表示觀測(cè)總數(shù)，Xi(i=1，…，n)表示第i個(gè)觀測(cè)的標(biāo)志物取值(不失一般性，假設(shè)標(biāo)志物值越高事件發(fā)生時(shí)間越短)。Ti(i=1，…，n)表示第i個(gè)觀測(cè)的結(jié)局事件理論發(fā)生時(shí)間。Zi=min(Ti，Ci)為觀察時(shí)間，其中Ci為刪失時(shí)間。δi=1(Ti≤Ci)為刪失指標(biāo)，其中1(·)為示性函數(shù)。

分析方法

1.問(wèn)題轉(zhuǎn)化

靈敏度Se(c，t)和特異度Sp(c，t)是依時(shí)ROC的決定因素。以CD定義為例，定義個(gè)體標(biāo)志物X取值大于閾值c為陽(yáng)性，則依時(shí)ROC下的靈敏度和特異度可分別表示為

Se(c，t)=P(X>c|T≤t)

(1)

Sp(c，t)=P(X≤c|T>t)

(2)

上式利用了貝葉斯定理進(jìn)行變換，目的在于將問(wèn)題轉(zhuǎn)化為給定X取值下的事件發(fā)生時(shí)間分布，以便更容易依據(jù)生存分析方法進(jìn)行估計(jì)(尤其是對(duì)刪失的處理)。可見(jiàn)對(duì)依時(shí)ROC指標(biāo)計(jì)算的關(guān)鍵在于估計(jì)X和T的聯(lián)合分布[8]，下文從非參、半?yún)⒎椒▋蓚€(gè)角度進(jìn)行總結(jié)。

2.非參方法

(1)Kaplan-Meier法

由Heagerty等人于2000年[5]提出。將式(1)和(2)改寫為生存函數(shù)形式

(3)

(4)

其中，S(t)為生存函數(shù)，S(t|X>c)是X>c時(shí)的條件生存函數(shù)，F(xiàn)X(c)是標(biāo)志物X的經(jīng)驗(yàn)分布函數(shù)，F(xiàn)X(c)=∑1(Xi≤c)/n。

對(duì)生存函數(shù)S(t)采用Kaplan-Meier法估計(jì)

(5)

Tn為觀察時(shí)間Zi值的集合。

Kaplan-Meier法簡(jiǎn)單易理解，但其獲得的靈敏度和特異度不具有單調(diào)性且不能確保落在[0，1]區(qū)間，違背上述指標(biāo)的基本特征。另外，生存函數(shù)的估計(jì)建立在刪失獨(dú)立于標(biāo)志物的條件下，該條件不滿足時(shí)會(huì)使得對(duì)結(jié)果的估計(jì)不穩(wěn)定。

(2)最近鄰估計(jì)法

考慮Kaplan-Meier法的上述問(wèn)題，Heagerty等人[5]同時(shí)提出了基于二元生存函數(shù)的最近鄰估計(jì)法(nearest neighborhood estimation，NNE)。假定刪失基于標(biāo)志物X條件獨(dú)立，從而采用S(t|X)而不是S(t)進(jìn)行靈敏度和特異度的估計(jì)。以NNE估計(jì)(X，T)的二元分布，條件生存函數(shù)的計(jì)算式為

(6)

于是，靈敏度和特異度為

Se(c，t)=P(X>c|T≤t)

(7)

(8)

(3)風(fēng)險(xiǎn)集遞歸算法

Chambless等人[9]于2006年基于生存分析中風(fēng)險(xiǎn)集的概念提出了一種類似于Kaplan-Meier法的遞歸算法。首先對(duì)觀察時(shí)間從小到大排序，令tk為排序后第k個(gè)觀測(cè)的觀察時(shí)間，tm為時(shí)刻t前的最后一次觀察到的事件的觀察時(shí)間。Blanche[10]等人在此基礎(chǔ)上給出了靈敏度和特異度更加直觀的表達(dá)。

(9)

(10)

其中，d(k)是在tk時(shí)發(fā)生事件個(gè)體的指標(biāo)，I(Xd(k)>c)=P(Xi>c|tk-1

與最近鄰估計(jì)法相比，風(fēng)險(xiǎn)集遞歸算法不涉及任何平滑參數(shù)直接得到生存函數(shù)，亦可保證靈敏度單調(diào)且范圍為[0，1]，但其特異度不單調(diào)且應(yīng)用于大型數(shù)據(jù)集時(shí)計(jì)算量大，故應(yīng)用并不廣泛。

(4)非條件/條件逆概率刪失加權(quán)法

逆概率加權(quán)的思想最早由Uno等人[11]于2007年引入依時(shí)ROC分析，經(jīng)Hung等人[12]于2010年發(fā)展，以非刪失概率的倒數(shù)(即逆概率)作為權(quán)重，故稱非條件逆概率刪失加權(quán)法(inverse probability censoring weighting，IPCW)。進(jìn)一步，Blanche等人[10]于2013年改用條件刪失逆概率作為權(quán)重，稱為條件IPCW。非條件IPCW的靈敏度為

(11)

其中，SC(Zi)是在觀察時(shí)間Zi上的刪失生存函數(shù)(即非刪失概率)。特異度為

(12)

將式(11)中的SC(Zi)更換為SC(Zi|Xi)，可得到條件IPCW的靈敏度估計(jì)

(13)

而此時(shí)的特異度也需做條件生存函數(shù)的轉(zhuǎn)化，為：

(14)

其中，SC(t|Xi)=P(Ci>t|Xi)。

刪失生存函數(shù)和條件生存函數(shù)均可采用KM估計(jì)，后者的優(yōu)勢(shì)在于其兼容不獨(dú)立于標(biāo)志物的刪失，即刪失與標(biāo)志物取值有關(guān)的情形。

(5)個(gè)體信息條件加權(quán)法

Li等人[13]于2016年提出一種容易理解且高效的加權(quán)算法。將研究對(duì)象在時(shí)刻t的狀態(tài)劃分為四組：對(duì)照(Zi>t)、病例(Zi≤t且δi=1)、此刻刪失(Zi=t且δi=0)、已經(jīng)刪失(Zi

Wi=P(Ti≤t|Zi，δi，Xi)

其中，S(t|X)=P(T>t|X)為條件生存函數(shù)，可用核加權(quán)KM估計(jì)。

靈敏度、特異度分別為：

(15)

(16)

該方法的優(yōu)勢(shì)在于：靈敏度和特異度關(guān)于閾值c單調(diào)；自動(dòng)解釋了刪失時(shí)間C與標(biāo)志物X之間可能存在的相關(guān)性；核加權(quán)KM估計(jì)對(duì)帶寬的選擇不敏感。

(6)加權(quán)平均秩法

基于ID定義和風(fēng)險(xiǎn)集個(gè)體得分的排序思想，Saha-Chaudhuri等人[14]于2013年提出加權(quán)平均秩法。該法的特點(diǎn)為不需計(jì)算靈敏度和特異度，而直接基于給定t時(shí)刻的局部一致估計(jì)得到依時(shí)ROC曲線下面積(area under curve，AUC)：

(17)

(18)

其中，Nt(hn)=(tj:|t-tj|

(19)

其中，Khn是標(biāo)準(zhǔn)化的核函數(shù)，滿足∑jKhn(t-tj)=1。式(19)是式(18)的平滑形式。該方法的假設(shè)為刪失時(shí)間C獨(dú)立于事件發(fā)生時(shí)間T。

(7)二元核密度估計(jì)法

CD定義下的靈敏度和特異度：

(20)

(21)

ID定義下的特異度與式(21)相同，而靈敏度：

(22)

3.半?yún)⒎椒?/p>

非參方法對(duì)函數(shù)形式假設(shè)條件較少，使用靈活，但逐點(diǎn)擬合過(guò)程計(jì)算量大。半?yún)?shù)法通過(guò)一定假設(shè)條件減少計(jì)算量，同時(shí)避免了參數(shù)法嚴(yán)格的條件限制[1]。

(1)非等比例Cox模型

ID定義與Cox比例風(fēng)險(xiǎn)模型中的風(fēng)險(xiǎn)函數(shù)直接對(duì)應(yīng)。由此Heagerty等人[6]2005年提出將比例風(fēng)險(xiǎn)回歸參數(shù)γ引入Cox模型，為了放寬Cox模型中等比例的假設(shè)條件，將固定比例γ擴(kuò)展為依時(shí)協(xié)變量系數(shù)γ(t)。

則靈敏度為：

(23)

特異度為：

(24)

(2)Cox平均生存概率法

一般而言閾值取值較大時(shí)對(duì)應(yīng)的樣本量較少，影響對(duì)條件生存函數(shù)S(t|Xi)估計(jì)穩(wěn)定性。Chambless等人[9]2006年針對(duì)該問(wèn)題提出了基于Cox模型“得分”的估計(jì)方法。其通過(guò)Cox比例風(fēng)險(xiǎn)模型估計(jì)風(fēng)險(xiǎn)因素的系數(shù)，進(jìn)而估計(jì)生存函數(shù)。其靈敏度和特異度分別為：

(25)

(26)

其中M=XTβ是生存函數(shù)的得分。該方法產(chǎn)生單調(diào)的敏感度和特異度，且刪失可以不獨(dú)立于標(biāo)志物取值。

(3)條件絕對(duì)風(fēng)險(xiǎn)函數(shù)法

Viallon等人[15]于2011年闡述了預(yù)測(cè)曲線與AUC的對(duì)應(yīng)關(guān)系，進(jìn)而提出一種基于條件絕對(duì)風(fēng)險(xiǎn)函數(shù)(conditional absolute risk function)直接估計(jì)依時(shí)AUC的方法，其通式為：

(27)

其中，F(xiàn)(t；Xi)=P(T≤t|X=Xi)為條件絕對(duì)風(fēng)險(xiǎn)函數(shù)，Xi表示排序后第i個(gè)觀測(cè)標(biāo)志物的值，可采用Cox比例風(fēng)險(xiǎn)模型、Aalen可加模型等方法進(jìn)行估計(jì)。

(4)分?jǐn)?shù)多項(xiàng)式估計(jì)法

基于ID定義，Shen等人[1]2015年提出一種采用分?jǐn)?shù)多項(xiàng)式直接估計(jì)AUC(t)的方法。令η(·)為連接函數(shù)，如logit函數(shù)，利用K階分?jǐn)?shù)多項(xiàng)式進(jìn)行建模

(28)

(5)貝葉斯半?yún)⒐烙?jì)法

Zhao等人[16]2016年采用線性相關(guān)的狄利克雷過(guò)程(the linear-dependent Dirichlet process，LDDP)，基于MCMC抽樣估計(jì)標(biāo)志物的分布函數(shù)和給定標(biāo)志物取值下的事件發(fā)生時(shí)間的條件分布。假設(shè)事件發(fā)生時(shí)間已經(jīng)轉(zhuǎn)換為對(duì)數(shù)形式。

CD和ID定義下靈敏度分別為：

(29)

(30)

其中，f(t|Xi)為條件密度函數(shù)，F(xiàn)(t|Xi)為條件累計(jì)分布函數(shù)，可通過(guò)LDDP混合模型計(jì)算得到。

CD、ID定義下的特異度均為：

(31)

其中，S(t|Xi)=1-F(t|Xi)為條件生存函數(shù)。

軟件實(shí)現(xiàn)

依時(shí)ROC分析可通過(guò)多種軟件實(shí)現(xiàn)。以R語(yǔ)言為例，依時(shí)ROC相關(guān)package見(jiàn)表1，讀者可在https：//cran.r-project.org檢索并查閱相應(yīng)幫助文檔。此外，DIVAT網(wǎng)站(http：//www.divat.fr/en/softwares)也提供了一系列用于特殊目的的依時(shí)ROC分析R package，如用于生物微陣列數(shù)據(jù)預(yù)后標(biāo)志物分析的ROC632。

表1 依時(shí)ROC分析方法R語(yǔ)言軟件包

趨勢(shì)和展望

近年來(lái)，依時(shí)ROC分析方法不斷與醫(yī)學(xué)研究數(shù)據(jù)和應(yīng)用需求結(jié)合，其方法學(xué)研究領(lǐng)域隨之不斷拓寬，包括但不限于以下發(fā)展趨勢(shì)：

1.重復(fù)測(cè)量數(shù)據(jù)縱向標(biāo)志物的依時(shí)ROC分析[17-19]；

2.考慮競(jìng)爭(zhēng)風(fēng)險(xiǎn)的生存分析評(píng)價(jià)[10，20-23]；

3.除常見(jiàn)的右刪失外，左刪失和區(qū)間刪失的依時(shí)ROC分析[23-24]；

4.標(biāo)志物存在缺失時(shí)性能的評(píng)價(jià)[25-26]；

5.特殊設(shè)計(jì)類型，如病例隊(duì)列研究、巢式病例對(duì)照研究下依時(shí)ROC分析[27-28]；

6.多標(biāo)志物的綜合評(píng)價(jià)[29-30]。

由于擴(kuò)展了經(jīng)典ROC曲線的時(shí)間維度，依時(shí)ROC在臨床應(yīng)用方面蘊(yùn)藏著巨大的價(jià)值[7]。陽(yáng)性預(yù)測(cè)值等指標(biāo)的估計(jì)方法[31]，及以最大化患者效用值為目標(biāo)的精準(zhǔn)分層方法[32]探討，一定程度上促進(jìn)了依時(shí)ROC方法向臨床應(yīng)用的轉(zhuǎn)化過(guò)程。然而，依時(shí)ROC分析方法的快速發(fā)展在一定程度上超過(guò)了醫(yī)學(xué)研究者的應(yīng)用速度，國(guó)內(nèi)研究中的應(yīng)用尚且較少[3，33-34]。如何針對(duì)所研究的問(wèn)題在眾多依時(shí)ROC分析方法中進(jìn)行選擇，仍需要開(kāi)展比較統(tǒng)計(jì)學(xué)研究。