本文原作者:尹迪,經授權后發布。

1 基本概念

1.1 生存數據

??生存數據就是關于某個體生存時間的數據。生存時間就是死亡時間減去出生時間。例如,以一個自然人的出生為“出生”,死亡為“死亡”。 那么,死亡時間減去出生時間,就是一個人的壽命,這是一個典型的生存數據。類似的例子,還可以舉出很多。所有這些數據都有一個共同的特點, 就是需要清晰定義的:出生和死亡 。如果用死亡時間減去出生時間,就產生了一個生存數據。因為死亡一定發生在出生的后面,因此,生存數據一定是正數。 因為,從理論上講,出生死亡時間都可能取任意數值,因此 生存數據一定是連續的正數。

??生存期不同于一般指標,他有二個特點:

• 1 有截尾數據(censored data)

??例如我們在疾病預測的實驗中,隨訪未能知道病人的確切生存時間,只知道病人的生存時間大于某時間。

(1)病人失訪或因其他原因而死亡---失訪 (2)到了研究的終止期病人尚未死亡---終訪

??例如,一個人的壽命。假設我關心1949年出生的人群的平均壽命。這群人可以被分成兩部分。一部分是已經離世了,所以他們的死亡時間是準確知道的。因此,他們的壽命是非常清晰的。 另一部分,是所有健在的人群,他們從1949年出生到現在,已經走過了將近70個春秋歲月,但是他們還活著!到2017年為止,他們已經生存了68年,但是他們最終的壽命是多少?我們是不知道的。 我們知道他們的壽命一定會比68大,數學上可以被記作68+。但是,到底“+”多少,不清楚。

??雖然截尾數據提供的信息是不完全的,但不能刪去,因為這不僅損失了資料,而且會造成偏性。

• 2 生存時間的特征一般不服從正態分布

??跟所有的數據分析一樣,要分析生存數據,首要問題是做描述性分析。如果生存數據沒有被截斷,那么所有常規的描述統計量,估計量都適用。例如:樣本均值,樣本方差等。 但是,如果生存數據存在大量的截斷數據,那么任何同均值相關的統計量就都沒法計算了。例如:樣本均值無法算,樣本方差涉及到因變量的平方的均值,因此它也沒法計算。

??真實的數據常常非常復雜,每個樣本的出生日期不同,死亡日期不同,截斷時間點不同。但是,不管這個數據如何復雜,其背后的基本原理是一樣的。 那就是:雖然樣本均值沒法估計,樣本方差沒法估計。但是,各種分位數卻在一個很大的范圍內可以被估計。如果這個范圍大到可以覆蓋中位數,那么從某種意義上講,我們也就把握了生存的平均狀況了。

??總結一下就是:對生存數據最基本的描述分析方法,不是過去常見的樣本均值,樣本方差等等,而是各種分位數。這些分位數也就構成了所謂的生存函數。生存函數就變成了對生存數據最基本的描述統計。

1.2 描述生存時間分布規律的函數

• 1 生存率(Survival Rate)

??又稱為生存概率或生存函數,它表示生存時間長于時間t的概率,用S(t) 表示:s(t)=P(T≥t)。以時間t為橫坐標,S(t)為縱坐標所作的曲線稱為生存率曲線,它是一條下降的曲線,下降的坡度越陡, 表示生存率越低或生存時間越短,其斜率表示死亡速率。

• 2 概率密度函數(Probability Density Function)

??其定義為:f(t)=lim (一個病人在區間(t,t+△t)內死亡概率/△t),它表示死亡速率的大小。如以t為橫坐,f(t)為縱坐標作出的曲線稱為密度曲線,由曲線上可看出不同時間的死亡速率及死亡高峰時間。 縱坐標越大,其死亡速率越高,如曲線呈現單調下降,則死亡速率越來越小,如呈現峰值,則為死亡高峰。

• 3 風險函數(Hazard Function)

??其定義為:h(t)=lim(在時間t生存的病人死于區間(t,△t)的概率/△t),由于計算h(t)時,用到了生存到時間t這一條件,故上式極限式中分子部分是一個條件概率。 可將h(t)稱為生存到時間t的病人在時間t的瞬時死亡率或條件死亡速率或年齡別死亡速率。當用t作橫坐標,h(t)為縱坐標所繪的曲線,如遞增,則表示條件死亡速率隨時間而增加,如平行于橫軸, 則表示沒有隨時間而加速(或減少)死亡的情況。

2 加速失效時間模型(AFT)

??在生存分析領域,加速失效時間模型(accelerated failure time model,AFT 模型)可以作為比例風險模型的替代模型。AFT模型將線性回歸模型的建模方法引人到生存分析的領域, 將生存時間的對數作為反應變量,研究多協變量與對數生存時間之間的回歸關系,在形式上,模型與一般的線性回歸模型相似。對回歸系數的解釋也與一般的線性回歸模型相似,較之Cox模型, AFT模型對分析結果的解釋更加簡單、直觀且易于理解,并且可以預測個體的生存時間。

??在spark ml中,實現了AFT 模型,這是一個用于檢查數據的參數生存回歸模型。它描述了生存時間對數的模型,因此它通常被稱為生存分析的對數線性模型。不同于為相同目的設計的比例風險模型(Proportional hazards model), AFT模型更容易并行化,因為每個實例獨立地貢獻于目標函數。

??給定給定協變量的值$x^{'}$,對于i = 1, …, n可能的右截尾的隨機生存時間$t_{i}$,AFT模型的似然函數如下:

$$L(beta,sigma)=prod_{i=1}^n[frac{1}{sigma}f_{0}(frac{log{t_{i}}-x^{'}beta}{sigma})]^{delta_{i}}S_{0}(frac{log{t_{i}}-x^{'}beta}{sigma})^{1-delta_{i}}$$

??其中,$delta_{i}$是指示器,它表示事件i是否發生了,即有無截尾。使$epsilon_{i}=frac{log{t_{i}}-x^{‘}beta}{sigma}$,則對數似然函數為以下形式:

$$iota(beta,sigma)=sum_{i=1}^{n}[-delta_{i}logsigma+delta_{i}log{f_{0}}(epsilon_{i})+(1-delta_{i})log{S_{0}(epsilon_{i})}]$$

??其中$S_{0}(epsilon_{i})$是基準生存函數,$f_{0}(epsilon_{i})$是對應的概率密度函數。

??最常用的AFT模型基于服從韋伯分布的生存時間,生存時間的韋伯分布對應于生存時間對數的極值分布,所以$S_{0}(epsilon)$函數為:

$$S_{0}(epsilon_{i})=exp(-e^{epsilon_{i}})$$

??$f_{0}(epsilon_{i})$函數為:

$$f_{0}(epsilon_{i})=e^{epsilon_{i}}exp(-e^{epsilon_{i}})$$

??生存時間服從韋伯分布的AFT模型的對數似然函數如下:

$$iota(beta,sigma)= -sum_{i=1}^n[delta_{i}logsigma-delta_{i}epsilon_{i}+e^{epsilon_{i}}]$$

??由于最小化對數似然函數的負數等于最大化后驗概率,所以我們要優化的損失函數為$-iota(beta,sigma)$。分別對$beta$和$logsigma$求導,得到:

$$frac{partial (-iota)}{partial beta}=sum_{1=1}^{n}[delta_{i}-e^{epsilon_{i}}]frac{x_{i}}{sigma}$$

$$frac{partial (-iota)}{partial (logsigma)}=sum_{i=1}^{n}[delta_{i}+(delta_{i}-e^{epsilon_{i}})epsilon_{i}]$$

??可以證明AFT模型是一個凸優化問題,即是說找到凸函數$-iota(beta,sigma)$的最小值取決于系數向量$beta$以及尺度參數的對數$logsigma$。 spark ml中使用L-BFGS作為優化算法。

注意:當使用無攔截(intercept)的連續非零列訓練AFTSurvivalRegressionModel時,Spark MLlib為連續非零列輸出零系數。這種處理與R中的生存函數survreg不同。

3 例子

 val dataList: List[(Double, Double, Double, Double)] = List(      (2, 51, 1, 1),      (2, 58, 1, 1),      (2, 55, 2, 1),      (2, 28, 22, 1),      (1, 21, 30, 0),      (1, 19, 28, 1),      (2, 25, 32, 1),      (2, 48, 11, 1),      (2, 47, 14, 1),      (2, 25, 36, 0),      (2, 31, 31, 0),      (1, 24, 33, 0),      (1, 25, 33, 0),      (2, 30, 37, 0),      (2, 33, 35, 0),      (1, 36, 25, 1),      (1, 30, 31, 0),      (1, 41, 22, 1),      (2, 43, 26, 1),      (2, 45, 24, 1),      (2, 35, 35, 0),      (1, 29, 34, 0),      (1, 35, 30, 0),      (1, 32, 35, 1),      (2, 36, 40, 1),      (1, 32, 39, 0))    val data = dataList.toDF("sex", "age", "label", "censor").orderBy("label")    val colArray = Array("sex", "age")    val assembler = new VectorAssembler().setInputCols(colArray).setOutputCol("features")    val vecDF: DataFrame = assembler.transform(data)    val aft = new AFTSurvivalRegression()    val model = aft.fit(vecDF)    // Print the coefficients, intercept and scale parameter for AFT survival regression    println(s"Coefficients: ${model.coefficients} Intercept: " +      s"${model.intercept} Scale: ${model.scale}")    val Array(coeff1, coeff2) = model.coefficients.toArray    val intercept: Double = model.intercept    val scale: Double = model.scale    val aftDF = model.transform(vecDF)    // 風險率h(t)    aftDF.selectExpr("sex", "age", "label", "censor",      "features", "round(prediction,2) as prediction",      s"round( exp( sex*$coeff1+age*$coeff2+$intercept ), 2) as h(t)").orderBy("label").show(100, false)

4 參考文獻

【1】Spark Doc

【2】回歸XY | 數據江湖:回歸五式之第五式(生存回歸)