依時ROC分析方法學綜述*

中國醫學科學院基礎醫學研究所/北京協和醫學院基礎學院統計學教研室(100005) 王子興 申郁冰 姜晶梅

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

本文從非參、半參角度總結依時ROC分析的方法學研究及其對應的軟件實現方式，以便增進讀者對依時ROC方法的了解并在臨床轉化中的應用提供參考。

基本定義

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

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

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

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

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

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

2.符號定義

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

分析方法

1.問題轉化

靈敏度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)

上式利用了貝葉斯定理進行變換，目的在于將問題轉化為給定X取值下的事件發生時間分布，以便更容易依據生存分析方法進行估計(尤其是對刪失的處理)。可見對依時ROC指標計算的關鍵在于估計X和T的聯合分布[8]，下文從非參、半參方法兩個角度進行總結。

2.非參方法

(1)Kaplan-Meier法

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

(3)

(4)

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

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

(5)

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

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

(2)最近鄰估計法

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

(6)

于是，靈敏度和特異度為

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

(7)

(8)

(3)風險集遞歸算法

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

(9)

(10)

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

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

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

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

(11)

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

(12)

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

(13)

而此時的特異度也需做條件生存函數的轉化，為：

(14)

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

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

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

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

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

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

靈敏度、特異度分別為：

(15)

(16)

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

(6)加權平均秩法

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

(17)

(18)

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

(19)

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

(7)二元核密度估計法

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

(20)

(21)

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

(22)

3.半參方法

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

(1)非等比例Cox模型

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

則靈敏度為：

(23)

特異度為：

(24)

(2)Cox平均生存概率法

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

(25)

(26)

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

(3)條件絕對風險函數法

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

(27)

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

(4)分數多項式估計法

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

(28)

(5)貝葉斯半參估計法

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

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

(29)

(30)

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

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

(31)

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

軟件實現

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

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

趨勢和展望

近年來，依時ROC分析方法不斷與醫學研究數據和應用需求結合，其方法學研究領域隨之不斷拓寬，包括但不限于以下發展趨勢：

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

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

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

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

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

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

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