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

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

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

基本定義

1.病例和對照的三種定義

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

圖1 依時ROC三種定義示意圖

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

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

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

2.符號定義

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

分析方法

1.問題轉(zhuǎn)化

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

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

(1)

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

(2)

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

2.非參方法

(1)Kaplan-Meier法

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

(3)

(4)

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

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

(5)

Tn為觀察時間Zi值的集合。

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

(2)最近鄰估計法

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

(6)

于是，靈敏度和特異度為

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

(7)

(8)

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

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

(9)

(10)

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

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

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

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

(11)

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

(12)

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

(13)

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

(14)

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

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

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

Li等人[13]于2016年提出一種容易理解且高效的加權(quán)算法。將研究對象在時刻t的狀態(tài)劃分為四組：對照(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估計。

靈敏度、特異度分別為：

(15)

(16)

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

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

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

(17)

(18)

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

(19)

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

(7)二元核密度估計法

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

(20)

(21)

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

(22)

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

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

(1)非等比例Cox模型

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

則靈敏度為：

(23)

特異度為：

(24)

(2)Cox平均生存概率法

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

(25)

(26)

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

(3)條件絕對風(fēng)險函數(shù)法

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

(27)

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

(4)分數(shù)多項式估計法

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

(28)

(5)貝葉斯半?yún)⒐烙嫹?/p>

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

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

(29)

(30)

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

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

(31)

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

軟件實現(xiàn)

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

表1 依時ROC分析方法R語言軟件包

趨勢和展望

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

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

2.考慮競爭風(fēng)險的生存分析評價[10，20-23]；

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

4.標志物存在缺失時性能的評價[25-26]；

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

6.多標志物的綜合評價[29-30]。

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