存活曲線估計 | 存活資料分析 | 生物統計學習筆記

存活分析（survival analysis）研究的是從明確起點到事件發生所經過的時間。事件可以是死亡，也可以是復發、痊癒、設備故障或其他事先定義的結果；因此「存活」在公式中代表事件尚未發生，不一定只指仍然存活。

母體存活函數與樣本存活曲線

令非負隨機變數 T 表示從研究起點到事件發生的時間。母體存活函數 S(t) 是事件時間超過 t 的機率，也就是個體在時間 t 之後仍未發生事件的機率。

母體存活函數

S(t)=P(T>t)=1-F(t)

沒有截尾時的樣本比例

\hat S(t)=\frac{\#\{T_i>t\}}{n}

沒有截尾時，Ŝ(t) 就是目前仍未發生事件的人數除以原始樣本數。但有截尾資料時，不能把離開研究的人算成事件，也不能在後續時間繼續把他放在分母；因此需要隨時間更新仍在風險中的人數。

右截尾（right censoring）表示在最後一次知道個體仍未發生事件的時間之後，真正事件時間未知。例如研究結束時仍未發生事件、失去追蹤或因其他原因退出，都可能形成右截尾。這筆資料仍告訴我們 T 大於已追蹤時間，所以不能直接刪除。

圖 A 以 calendar time 顯示不同受試者在不同日期進入研究；直接比較線段終點會混合入組日期與實際追蹤長度。圖 B 把每個人的入組時間對齊為 t=0，橫軸才是 survival time。實心圓代表觀察到事件，空心圓與短線代表右截尾。

將不同事件時間由小到大記為 t₁,t₂,…。在第 j 個事件時間發生前，令 nⱼ 為仍在風險集的人數，dⱼ 為該時間發生的事件數。從該時間點存活過去的條件比例是 (nⱼ−dⱼ)/nⱼ。

第 j 個事件時間的條件存活比例

P(T>t_j\mid T\ge t_j)=\frac{n_j-d_j}{n_j}

Kaplan–Meier product-limit estimator

\hat S(t)=\prod_{t_j\le t}\left(\frac{n_j-d_j}{n_j}\right)

存活超過後一個時間點，必須先存活超過前面的所有事件時間，所以各段使用乘法而不是加法。每發生一次事件，就在前一階的累積存活率上乘上新的條件存活比例，形成向下的階梯曲線。

把圖中的追蹤時間由小到大排列。時間後面的「+」表示右截尾；截尾列不建立新的存活乘數，因此累積存活率維持不變。下表使用與重製圖相同的資料。

t=2

\hat S(2)=\frac9{10}=0.900

t=6

\hat S(6)=\frac9{10}\frac8{9}=0.800

t=7：同時發生兩個事件

\hat S(7)=0.800\times\frac{8-2}{8}=0.600

t=8

\hat S(8)=0.600\times\frac45=0.480

t=9

\hat S(9)=0.480\times\frac34=0.360

t=12

\hat S(12)=0.360\times\frac12=0.180

t=7 發生 A、C 兩個事件後，風險集由 8 人減為 6 人；I 又在 7+ 被截尾，所以進入 t=8 前只剩 5 人。E 在 11+ 被截尾後，t=12 的風險集只剩 B、D 兩人。這正是不能始終用原始 10 人作分母的原因。

Ŝ(t) 是由樣本估計的階梯曲線，不同樣本會得到不同結果。Greenwood 公式把每個事件時間的變異貢獻 dⱼ/[nⱼ(nⱼ−dⱼ)] 累加，再乘上目前累積存活率的平方，得到 Ŝ(t) 的估計變異數。

Greenwood variance

\widehat{\operatorname{Var}}\!\left[\hat S(t)\right]=\hat S(t)^2\sum_{t_j\le t}\frac{d_j}{n_j(n_j-d_j)}

Standard error

SE\!\left[\hat S(t)\right]=\hat S(t)\sqrt{\sum_{t_j\le t}\frac{d_j}{n_j(n_j-d_j)}}

原始表格中的 standard deviation 更精確地說是 standard error，因為它描述的是估計曲線 Ŝ(t) 在重複抽樣下的不確定性，而不是個別受試者存活時間的標準差。

一般 Wald 型信賴區間

\hat S(t)\pm z_{1-\alpha/2}SE\!\left[\hat S(t)\right]

95% 信賴區間

\hat S(t)\pm1.96\,SE\!\left[\hat S(t)\right]

令 g(t)=ln[−ln Ŝ(t)]。只要 0<Ŝ(t)<1，g(t) 可以落在整條實數線上；先在 g(t) 的尺度建立對稱區間，再轉回存活機率，就能自然把上下界保留在 0 與 1 之間。

轉換後的標準誤

SE\!\left\{\ln[-\ln\hat S(t)]\right\}=\sqrt{\frac{1}{[\ln\hat S(t)]^2}\sum_{t_j\le t}\frac{d_j}{n_j(n_j-d_j)}}

先在 log–log 尺度建立區間

g(t)\pm z_{1-\alpha/2}SE\{g(t)\}

轉回存活率尺度

\hat S(t)^{\exp[z_{1-\alpha/2}SE\{g(t)\}]}<S(t)<\hat S(t)^{\exp[-z_{1-\alpha/2}SE\{g(t)\}]}

存活函數回答「到時間 t 仍未發生事件的機率」；危險函數 h(t) 則回答「已經存活到 t 的個體，在接下來極短時間內發生事件的瞬時速率」。分母必須以目前仍在風險中的個體為條件，而不是以最初全體為分母。

Hazard function

h(t)=\lim_{\Delta t\to0}\frac{P(t\le T<t+\Delta t\mid T\ge t)}{\Delta t}

與機率密度及存活函數的關係

h(t)=\frac{f(t)}{S(t)}

由 F(t)=1−S(t) 得到

f(t)=F'(t)=-S'(t)\quad\Longrightarrow\quad h(t)=-\frac{S'(t)}{S(t)}

累積危險函數

H(t)=\int_0^t h(u)\,du=-\ln S(t)\quad\Longleftrightarrow\quad S(t)=e^{-H(t)}

f(t) 描述事件時間落在 t 附近、相對於原始母體的密度；除以 S(t) 後，才轉成以已經存活到 t 的風險集為基準的瞬時事件率。這也說明了累積危險 H(t)=−ln S(t)，並連回前面 log–log 信賴區間使用的轉換。