動態事件時間數據的多任務Logistic生存預測方法

阮燦華，林甲祥

（福建農林大學計算機與信息學院，福州350002）

（?通信作者電子郵箱ruancanhua@fafu.edu.cn）

0 引言

生存分析（survival analysis）是一種估計事件時間（timeto-event）的統計方法，已被廣泛應用于臨床醫學、衛生保健、生物信息學、人口統計學、機械制造、經濟學等領域［1］。在生存分析問題中，“生存”不僅可以指人或動物存活的狀態（相對于死亡事件），還可以表示某藥物的有效作用（相對于失效事件）、患者病情穩定的狀態（相對于疾病復發或病情惡化事件），以及儀器設備的正常工作狀態（相對于故障事件）等。因此，“生存時間”也是廣義的，不單單是指通常意義下生物體的存活時間，而是泛指研究者所關心的某現象（如患者體征穩定）的持續時間。生存分析的核心問題是對事件時間進行估計，即預測某研究對象的生存時間或特定時間段內某事件未發生的概率。生存分析是研究特定人群壽命、各類慢性疾病、藥物療效以及提高醫療衛生服務的一項重要課題。例如，臨床醫學研究需要預估冠心病患者出院后下一次可能復發的時間；流行病學和生物學研究則需要檢測樣本從接觸感染到發病所需要的時間；機械制造過程中人們常常對設備的可靠性進行測試，并估算故障發生時間（即設備使用壽命）。

傳統生存分析方法按照生存時間分布假設及模型參數可以分為無參、有參和半參三種類型。無參模型主要包括Kaplan-Meier（KM）估計［2］、Nelson-Aalen 估計［3］以及生命表法（life table）；有參模型通常假設事件時間服從某種特殊分布時，并采用正態、韋伯（Weibull）、log-Logistic等分布函數估計生存概率。然而，大多數情況下，關于事件時間的先驗知識相對較少，因此應用局限性較大。半參（semi-parametric）模型的主要代表是比例風險（proportional-hazard）Cox回歸模型［4］，該模型具有簡單易用的特點。近年來，機器學習方法也被廣泛用于生存數據處理以及生存模型的參數學習：Leblanc等［5］提出了一種針對刪失數據的分類回歸樹，即生存樹（survival tree）；Khan等［6］提出了一種針對刪失數據的支持向量回歸（Support Vector Regression，SVR）模型；楊銘等［7］針對復雜中醫腫瘤臨床數據提出了一種采用生存網絡結構模型；潘浩等［8］提出了一種基于深度學習的Cox-Lasso生存分析方法用于肺癌細胞自動檢測；Fisher等［9］針對縱向數據（longitudinal data）提出了一種動態Cox模型；Cortese等［10］提出了一種可以處理內在（internal）風險因子變量的多狀態對比風險模型；Wang［11］提出了一種針對動態風險因子的數值插入法以及對未知鏈接函數進行估計的動態比例風險模型。

研究事件發生規律常常需要對樣本進行長期追蹤，因此本文需要處理的事件時間數據包含大量隨時間變化的風險因子變量（如哮喘病患者在不同時間點的血液含氧量和二氧化碳分壓觀測值）。然而，現有的生存分析方法大多無法有效處理這類復雜的動態數據，這些方法通常采用簡單的單點取樣［12］（如當前值抽樣和平均值抽樣），在估計某時刻的生存概率時只考慮當前記錄時間點的風險因子變量值而忽略歷史數據的作用，這在許多實際應用中并不合理，例如哮喘患者在服用β2激動劑之后，本文會觀測到生理指標一秒鐘用力呼氣容積（Forced Expiratory Volume in one second，FEV1）在很長一段時間內都會受該藥劑影響，因此某一時刻的FEV1觀測值并不能刻畫病人的病情發展以及變量FEV1與生存概率之間的關系?？傮w而言，如何有效分析復雜的動態事件時間數據、構建事件時間預測模型是一個亟待解決的難點和熱點問題，因此，本文的工作重點是探究刻畫風險因子變量觀測序列的數據模型以及生存機器學習方法。

事件時間數據作為一種特殊的觀測型數據，包含大量刪失數據，可以用于事件預測模型學習的事件時間信息多有缺失。多任務學習方法能夠為每個訓練任務提供更多的信息，從而提高模型精度并降低預測錯誤率，例如，多任務學習預測病情發展、HIV（Human Immunodeficiency Virus）治療、基因數據分析等［13-14］。鑒于此，本文提出了一種多任務Logistic生存學習模型（Multi-task Logistic Survival Learning，MLSL），該模型將生存預測任務轉化為一系列不同時間點的Logistic生存回歸學習與分類問題，即多任務學習與預測。針對風險因子變量的動態觀測值，本文將構建一個新的全數據生存函數，利用所有可觀測值估計累計事件風險，提高數據利用率；另外，本文的新方法利用Logistic回歸模型分別估計事件樣本和刪失樣本在不同時間點的事件概率，并根據各時間點的概率對事件時間進行估計。為了學習MLSL模型的參數，本文在優化目標函數中對模型參數進行了正則化和平滑處理，避免模型過擬合。最后，本文在多個臨床事件時間數據上開展對比實驗，并驗證模型的有效性和實用性。

1 生存預測

本章介紹事件時間數據，并描述動態事件時間數據特性以及Logistic生存函數和概率模型。

1.1 事件時間數據

不同于一般類型的數據，事件時間數據是一種特殊的觀測型數據，常常包含右刪失（right-censoring）、左刪失、區間刪失等類型的事件刪失樣本。事件刪失的主要原因有以下兩種：

失訪 某些與研究無關的原因導致樣本無法訪問，如肺癌患者在隨訪期間死于心肌梗塞、自殺或車禍。

追蹤終止 追蹤在事件發生之前終止。

如圖1所示的7個樣本中，A、B、G的初始觀測時間未知，即事件左刪失；樣本B、C、E中途退出實驗，G在追蹤終止時仍未有事件發生，這4個樣本都屬于事件右刪失類型。本文主要針對最常見的右刪失類型數據開展模型研究和應用研究。

圖1 事件時間數據示例Fig.1 An exampleof time-to-event data

1.2 動態事件時間數據

樣本i在t時刻的D個風險因子變量的觀測值記為x（i t）=（xi（1t），xi（2t），…，xiD（t））。對于給定數據集，假設所有非刪失樣本的事件時間點為τ1＜τ2＜…＜τK。給定風險因子變量觀測時間點集合為φ，動態事件時間數據X i=（x（i u1），x（i u2），…，xi(u|φ|)），其中uj∈φ是風險因子變量觀測時間點。特別地，本文用φ（t）表示在時間t之前的所有觀測時間點。生存狀態zi為表示事件是否發生的二元變量，即zi∈｛0，1｝：當i的事件已發生則取值為1，否則為0；標簽時間變量yi∈R+在zi=1時為生存時間Ti，即yi=Ti，而當zi=0，則為刪失事件時間Ci，即yi=Ci。本文用E1={i|zi=1}表示追蹤期內有事件發生的所有樣本。對于刪失事件樣本，本文無法在追蹤期內觀測到事件是否發生，也無法得知其在追蹤期外的生存狀態，但至少本文已知在丟失其可觀測的最后時間點之前仍是“存活”的。換言之，右刪失樣本的生存時間大于或等于刪失時間，即Ti≥Ci，?i∈E0，集合E0=｛i|zi=0｝包含所有追蹤期內沒有事件發生的樣本。

1.3 Logistic生存函數

生存預測的核心任務是預測事件時間，即預測事件在某時間點發生的概率，其概率密度函數為F(t)=Pr(T≤t)。生存函數S（t）表示t時刻的生存概率，即事件不會早于t發生的概率，因此S(t)=Pr(T＞t)=1-F(t)。實際上，任意形式的分布（如指數分布、Weibull分布、正態分布），只要滿足時間變量t∈［0，∞）都可以作為生存分布假設［15］。例如當且僅當lnT=η+σU時（U∈（-∞，∞）為隨機變量），T服從log-Logistic分布，生存概率則為：

從式（1）可以看出，生存函數與風險因子變量無關，因此傳統的log-Logistic有參模型無法估計風險因子變量對事件發生的影響。本文改寫η為風險因子變量的線性函數，并構造如下生存函數：

這里，參數w刻畫了風險因子自變量在t時刻對事件發生的瞬時作用。生存函數S表示觀察對象隨訪到t時刻的累積生存率，該生存概率只由當前觀測值決定，而風險因子觀測值的變化對事件時間的影響作用被忽略。

2 多任務Logistic事件時間預測

本章給出新的風險累積Logistic生存函數和多任務學習方法、回歸模型系數的選擇策略，以及事件時間估計方法。

2.1 風險累積Logistic生存函數

為了估計風險因子變量在整個觀測期內對生存概率的累積作用，本文構建以下Logistic生存函數：

其中Ri（u，W）是i在觀測時間點u的累積事件風險，計算如下：

指數函數η（u，t）為u時刻的事件風險到t時刻的衰減率，取值范圍為（0，1］，即時間越久，風險衰減越多。換言之，累積生存函數可以充分利用觀測數據，更符合臨床應用中的實際情況，例如病人用藥之后的很長一段時間內，生理指標都會受到藥效影響，且該影響會隨著時間變小。

2.2 多任務生存學習

實際上，風險因子變量與事件發生概率之間的關系是隨時間而變化的，因此本文設定一個新的回歸系數矩陣W=（w0，w1，w2，…，w K）用于描述風險因子變量對生存概率的潛在影響力：w k是在tk時刻的D維回歸系數向量。樣本i的生存狀態在觀測期內的變化可以表示成一個二元序列，zi（τ2），…，zi（τK）），當t≤Ti時，zi（t）取值為0，t＞Ti時則取值為1。為了求取最優的回歸系數W*，本文將優化關于W的似然函數。該似然估計實際上是在給定參數時對z i概率的估計。給定i∈E1，關于W的似然函數為：

最大化似然函數可以轉化為最小化似然函數的負對數，由此，本文的學習任務為：

2.3 正則化系數選擇

正則化系數λ可以在驗證集上采用交叉驗證方法進行優化。為了避免優化過程中目標函數出現常數項，本文對訓練集風險因子變量矩陣X和驗證集風險因子變量矩陣U的各列，以及訓練集事件輸出列向量Y和驗證集事件輸出列向量Q都進行標準化處理。最優的λ值能夠使得模型在驗證集上的預測錯誤率最低，因此正則化系數優化的核心是解決下述最小化問題：

為了計算方便，上述優化問題的解可以近似估計為ridge回歸模型的正則化系數最優解，因此

目標函數O關于λ的偏導數為：

根據推論1以及矩陣計算，本文可以得到：

將式（11）和式（9）代入式（10），并令J=XTX+λI，則本文有：

式（12）包含大量的矩陣轉置計算，難度大且效率低，因此，本文用XTX的分解特征值來改寫式（12）。假設矩陣V的每列為特征向量，γ為特征值向量，那么XTX=Vdiag（γ）VT。G=XTX+λI的特征向量仍然為V，特征值則為γ+λI。因此，本文有G=Vdiag（γ+λI）VT。假設P=diag（（γ+λI）-1），P-1=VPVT，則

推論1 假設M是一個基于實參κ的方陣，且各分量函數可導，M（κ）對于所有的T都是可逆的，那么d

證明假設mij（κ）是M的分量函數，mjk（κ）是M-1（κ）的分量函數。給定矩陣M的階n以及克羅內克（Kronecker）函數，對于任意κ，本文有：

其中Ξ=｛τ1，τ2，…，τK｝。損失函數Δ（G，T）為預測事件時間G和真實事件時間T之間的誤差，如絕對誤差|G-T|或對數誤差|lnG-lnT|。

2.5 模型分析

本文的多任務回歸學習方法將生存函數（曲線）的擬合問題轉化為各時間點的生存狀態預測問題。這些預測任務是緊密相關的，各時間點的預測結果滿足生存概率的單調性。

與參數模型相比，本文的新模型無需生存時間分布假設，而且可以給出各時間點的事件風險率以及風險因子變量對風險率的影響，具有一定的普遍適用性。Cox模型采用比例風險假設，并估計生存概率為：

其中：基線風險H（0t）為不考慮風險因子變量時的累積風險，通常由Breslow估計法［16］確定；基線概率SCo（xt）=S（t（|0，0，…，0）），表示在t時刻的基準死亡比例（即發病密度或死亡密度）?；貧w系數β描述了生存概率與風險因子變量之間的關系，該關系不隨時間變化而變化。靜態回歸系數簡單，但卻不如本文的動態回歸參數更符合實際情況，這也是動態Cox（Time-dependent Cox，TCox）模型［9］采用動態系數βk的主要原因。類似地，Cox的半比例風險（Semi-Proportional Hazards，SPH）［17］模型也采用動態回歸系數。與SPH模型的全數據學習方法類似，新模型充分利用刪失數據優化參數。本質上，新模型和基于Cox的模型采用的都是風險乘積方式。

新模型中的全數據似然函數不同于Cox模型的最大化如下局部似然估計（Maximum Partial Likelihood Estimate，MPLE）［18］函數：

2.4 事件時間預測

模型學習在訓練樣本上優化生存函數參數并得到最優值W*。因此，對給定風險因子變量值Xtest的測試樣本，本文用Logistic回歸函數估計生存概率為S（t|Xtes）t=（1+R（i u，W*））-1。鑒于生存函數只能估計生存概率，且目前尚沒有具體的生存時間估計方法，本文對樣本i的事件時間估計如下：

其中χ（yi）是在時間yi前所有未有事件發生的樣本。另外，新模型的目標函數沒有采用傳統的Elastic-Net［19］懲罰因子

3 實驗

本章首先說明數據集及預處理工作，接著給出對比模型和性能評價標準，最后介紹實驗過程并根據實驗結果對模型的有效性和效率進行分析。

3.1 實驗數據及預處理

本實驗將在如表1所示的10個數據集上開展，其中，Lung、Prostate和Colorectal三個癌癥數據集是由NIH CDAS（https：//cdas.cancer.gov）提供；其他7個數據集從R database獲得，且本文利用tmerge函數［20］為每個數據集生成動態風險因子變量值。為了減少觀測值差異對相似度度量的影響，本文將數值型數據集的觀測值都進行標準化處理。對于類屬性變量，本文利用虛擬變量法（dummy variable）將其轉化為0-1二元變量，這樣所有風險因子變量便可以直接用于回歸模型的計算。虛擬變量法的核心是將有l(l≥2)個觀測值的分類型風險因子變量用其中任意l-1個觀測值所構成的0-1二元虛擬變量所表示。例如，生理指標BMI（Body Mass Index）一般有“體重過輕”“正?！薄俺亍焙汀胺逝帧彼膫€類屬性變量值，本文任選三個變量值作為虛擬變量來描述變量BMI。假設（體重過輕、正常、超重）為一組虛擬變量，每個虛擬變量取值1（“是”）或0（“否”），（0，0，0）則表示“肥胖”。本文都從每組數據中抽取10%的樣本作為驗證集用于選擇MLSL的正則化系數。

表1 事件時間數據集統計信息Tab.1 Statisticsof time-to-event datasets

3.2 對比模型

本文將MLSL與以下6個具有代表性的生存預測模型作對比：

1）TCox，即動態 Cox模型［9］，由 R 提供的 survival宏包實現。

2）CoxEN，即彈性網絡（Elastic-Net）規范化Cox模型［19］，由R宏包survival和glmnet實現。

3）RSF，即隨機生存森林（Random Survival Forest）［21］，可以用randomForestSRC包實現，模型參數為R函數默認值。

4）Cox-TRACE，即基于Cox的跡準則多任務生存分析模型，由文獻［22］提供的開源代碼實現。

5）SPH，即半比例風險模型，以Cox為基礎用于處理序列型特征變量，本文對比實驗采用與文獻［17］相同的參數設置。

6）AFT，即疊加失效時間（Accelerated Failure Time）模型［23］，假設事件時間滿足 Weibull分布，由R函數WeibullReg實現。

3.3 評價指標

為了評價各模型的生存概率預測能力，本文采用常用的AUC（Area Under the Curve）和C-index［24］作為評價指標。在隨訪學習結束時，預測事件是否發生是一個二元分類問題，因此本文采用如下的AUC指標檢測模型的預測（二分類）能力：

其中：I是0-1示性函數，t0為隨訪時間。

C-index指標是AUC指標的泛化形式，用于評價模型概率預測的能力。若數據集共有n個可以按照真實事件時間排序的樣本組合，則：

其中：i∈E1，yi≤yj表示所有可排序的樣本組合i、j。C-index統計在各觀測時間點yi所有滿足預測事件概率的排序與真實事件時間排序一致的樣本對i、j的個數，并考慮了生存并列（ties）樣本的存在［25］。AUC和C-index指標的取值范圍通常是0.5（完全隨機預測）到1.0（最佳預測）。

為了檢測模型預測事件時間的準確率，本文對比預測時間G和真實事件時間T之間的誤差，并定義如下相對錯誤率（Relative Error rate，RE）：

顯然，RE值越小，模型預測的事件時間越準確。

3.4 結果與分析

3.4.1 生存預測結果對比

每個模型在10個數據集上各進行100次10折交叉驗證，所得結果的平均值為最終結果。表2的每個單元格列出了平均AUC和平均C-index值，同時表3給出了對應的MLSL模型正則化系數。

表2 各模型關于AUC和C-index指標的預測結果比較Tab.2 Comparison of prediction resultsof each model in termsof AUCand C-index

表3 MLSL正則化系數Tab.3 Regularized coefficientsfor MLSL

表2結果表明不論是對于刪失比例很高的PBC、Nwtco、Veteran和Lung四個數據集，還是在樣本數和維度都較大的Prostate和Colorectal數據集上，MLSL均獲得了較高的預測精度，由此可見本文提出的MLSL模型對各類復雜數據具有較強的適用性。盡管對比模型能夠在個別數據集上取得較好的結果，如SPH在CGD和Lung數據集上，Cox-TRACE在Standford2數據集上達到或超越了MLSL的預測精度，但是整體表現不如MLSL，其主要原因就是這兩個模型是基于Cox的等比例風險假設。同樣地，TCox和CoxEN魯棒性也較差。參數模型AFT在所有數據集上都表現較差，原因在于Weibull分布假設并不符合大多數數據的真實情況。值得注意的是，多數情況下，模型在同一數據集上獲得的C-index精度往往比AUC略低，這是因為C-index作為泛化的AUC，在計算過程中需要對全部樣本生存時間進行精確排序，難度較大。

表4給出了模型在測試集上的事件時間預測結果的相對錯誤率。除了在FLchain數據集上相對錯誤率比SPH的稍大以外，MLSL關于事件時間的預測結果都是最好的，這主要得益于全數據學習和多任務學習對生存概率的估計和生存函數參數的優化。另外，生存函數本身也對事件時間的預測起決定性作用，本文的Logistic生存函數只需要在參數得到優化的情況下就可以較精確地擬合數據集中的樣本生存概率分布。然而，基于Cox的方法雖然在參數學習過程不需要估計基線風險，但是在預測生存概率和事件時間時則需要用到某些參數估計方法（如Breslow估計）［16］，從而影響了預測準確率。

3.4.2 生存概率曲線對比

圖2繪出了一位67歲的男性肺癌（Lung）患者和一位82歲女性直腸癌（Colorectal）患者在一年內（52周）的生存概率曲線。肺癌患者和直腸癌患者的生存時間分別為13周和10周，因此從理論意義上說，從第13周開始，肺癌患者的生存概率為0，即S（t＞13）=0，而直腸癌患者則從第10周開始生存概率應當為0，即S（t＞10）=0。從圖2（a）可以看出MLSL在第13周時所預測的生存概率最低，且早在第8周，MLSL預測該肺癌患者的生存概率低于50%，遠比其他預測模型所得到的生存概率低。在圖2（b）中，MLSL預測的結腸直腸癌患者生存概率在第10周之后是最低的。由此可見，MLSL模型具有一定的預警能力，能夠盡早發出警告，這在臨床決策具有相當重要的意義。另外，目標函數中的參數平滑處理使得MLSL生存曲線要比其他曲線略微平滑。

表4 事件時間預測的相對錯誤率比較Tab.4 REcomparison of predicted time-to-events

圖2 各模型預測的生存概率曲線對比Fig.2 Comparison of predicted survival curves by different models

3.4.3 模型效率對比

對比實驗中預測模型的測試效率都較高，且差異不明顯，因此本文只給出各模型的訓練效率。圖3列出了各模型在10個數據集上的訓練時間（單位：s）。由圖可見，MLSL在大數據集上的效率高于同樣采用動態系數的TCox模型，這主要是由于新模型的似然估計對每個樣本的生存概率估計是相互獨立的，目標函數收斂速度比偏似然估計快。CoxEN需要估計兩個模型參數，當風險因子變量較多時，參數優化消耗時間較多，因此處理維度較高的數據時效率較低。RSF模型中，樹節點的分裂過程需要評估log-rank結果，訓練時間受樣本數量影響較大，因此處理不同規模數據的效率差異明顯。SPH優化具有瞬時風險約束條件的目標函數，模型復雜度相對較高，算法迭代次數增加，因而其運行時間也有所增加。

為了進一步分析MLSL的泛化能力，本文檢測了MLSL的效率與樣本數量N和數據維度D之間的關系。以Prostate數據為例，本文重新合成該數據集，分別調整N和D并測試MSLS的效率：固定數據維度D=30，并不斷增加樣本數量，MLSL的運行時間與樣本數呈線性關系，如圖4（a）所示；當固定樣本數量N=3500，并不斷增加風險因子個數，MLSL的運行時間與維度D也呈現線性關系，如圖4（b）所示。以上結果表明本文的MLSL具有很強的泛化能力，適用于不同規模的數據。

圖3 MLSL的訓練時間與數據規模和數據維度之間的關系Fig.3 Relationship between training time and data size/dimensionality of MLSL

4 結語

本文提出了一種針對動態事件時間數據的多任務Logistic生存學習模型MLSL。新模型通過全數據學習方式估計事件樣本和刪失樣本的生存概率，并構建了具有參數正則化約束和平滑約束的優化目標函數。MLSL將生存預測轉化為一系列在不同時間點的Logistic回歸學習任務，無需未知數據生存分布假設，同時增強了模型對刪失數據的預測能力。針對動態數據的特性，MLSL刻畫了風險因子序列特征，并估計事件累積風險。與傳統方式相比，新模型能夠挖掘出風險因子變量與生存概率之間的潛在關系。在實際臨床數據上的對比實驗結果表明，MLSL可以有效地預測動態事件時間數據，其預測精度優于目前經典的生存分析方法，并且大幅度提升了預測結果的穩定性。