參考資料:
1. Competing Risks(Germ´an Rodr´ıguez)
2. Competing Risks - What, Why, When and How? (Sally R. Hinchliffe)
1. 簡介與例子引入
competing risks指造成failure的原因有多種可能。
以放置節育環為例，這里的failure指的應該是節育環沒能在原定期限前一直發揮好作用(即避孕作用)。導致failure的risk包括意外懷孕、節育環排出、醫學原因移除、個人原因移除等；由不同原因造成的failure，視為不同類的failure。
對於這類問題，研究的興趣點主要在於以下三項:
1) 研究變量向量$x$和特定failure的risk(指造成failure的為特定原因的風險risk)之間的關系；
2) 分析在某類failure高風險的人群是否也存在另一類failure的高風險；
3) 估計某類failure在排除其他原因的情況下的risk。
如果有不同類failure相互獨立的假設，問題3)可以被回答被計算(獨立的話，條件概率和自己本身的概率是一樣的)，而問題2)會因此得到否定的回答。
2. 標記說明和基本公式
$T$: Survival Time(生存時間)
$J$: The type of failure(failure的類型，即因為何種原因造成的failure).$j\in \{1,2,...,M\}$
$x$: A vector of covariate(協變量的向量)
總體風險比率(overall hazard rate)仍像一般情形的定義:
\[\lambda(t,x)=\lim_{\Delta t \to 0}\frac{P\{t\leq T < t+\Delta t|T\geq t, x\}}{\Delta t}\]
特定原因導致failure的瞬時風險(cause-specific hazard rate)則為:
\[\lambda_j(t,x)=\lim_{\Delta t \to 0}\frac{P\{t\leq T < t+\Delta t, J=j|T\geq t, x\}}{\Delta t}\]
這里假設failure只能由其中一種原因導致。如果存在兩種以上failure同時發生(同時有兩種最開始定義的原因導致了failure)，把這種情況定義成新的一類failure(兩種原因同時出現時，定為新的原因，和原來的兩種獨立造成failure的情況區分開來)，使假設依然成立。
這樣子，就有
\[\lambda(t,x)=\sum_{j=1}^{m} \lambda_j(t,x)\]
根據上式及生存函數的定義得總體生存函數(overall survival function)為
\[S(t,x) = e^{-\Lambda(t,x)} = e^{-\int_0^t \lambda(u,x)du}\]
(當$x$為time-varying變量時，函數形式會有變化)
類似地，定義
\[S_j(t,x) = e^{-\Lambda_j(t,x)},\]
\[\Lambda_j(t,x)=e^{-\int_0^t \lambda_j(u,x)du}\]
$*$當failure類型數$M>1$時，一般來說，$S_j(t,x)$不能被理解成生存函數survival function，但可以推出$S_j(t,x)$和$S(t,x)$之間的關系:
\begin{eqnarray*}
S(t,x) & = & e^{-\int_0^t \lambda(u,x)du} \\
& = & e^{-\int_0^t \sum_j \lambda_j(u,x)du} \\
& = & e^{-\sum_j \int_0^t \lambda_j(u,x)du} \\
& = & \prod_j e^{-\int_0^j \lambda_j(u,x)du} \\
& = & \prod_j S_j(t,x)
\end{eqnarray*}
同樣地，定義cause-specific的t時刻密度函數:
\[f_j(t,x)= \lim_{\Delta t \to 0} \frac{P\{t\leq T \leq t+\Delta t, J=j|x\}}{\Delta t}\]
\[f_j(t,x)= \lambda_j(t,x)S(t,x).\]
這是衡量個體在t時刻因原因$j$死亡的非條件概率風險，並且有
\[f(t,x) = \sum_{j=1}^M f_j(t,x)\]
3. 估計(one sample:同質樣本，無協變量，非參情形)
3.1 Kaplan-Meier
令$t_{j1} < t_{j2} < \cdots < t_{jk_j}$，表示的是第$j$類failure的failure time。
$n_{ji}$表示在$t_{ji}$因原因$j$死亡的個體數。
用一般KM估計量的推導方式可以得到KM估計量
\[\hat{S_j}(t,x)=\prod_{i:t_{ji}} (1-\frac{d_{ji}}{n_{ji}})\]
可以看到，對於某個$j$類failure，其他類failure的個體的處理和出現右刪失的個體的處理結果是一樣的。
如果不同類failure之間沒有結(ties)，則
\[\hat{S}(t) = \prod_{j=1}^{M} \hat{S_j}(t),\]
$S(t)$的估計量是cause-specific survivor-like functions(指$\hat{S_j(t)}$)的估計量的乘積。
3.2 Nelson-Aalen
cause-specific累積風險函數的估計:
\[\hat{\Lambda_j}(t)=\sum_{i:t_{ji}<t} \frac{d_{ji}}{n_{ji}}\]
對應着
\[\hat{\lambda_j}(t)=\frac{d_{ji}}{n_{ji}},t=t_{ji}.\]
基於K-M估計的另一種估計:
\[\hat{\Lambda_j}(t)=-\log \hat{S_j}(t)\]
同理也可以從$\hat{\Lambda_j}(t)$的估計得到$\hat{S_j}(t)$的另一種形式的估計。
4. 估計:引入回歸模型
n個observation:$(t_i, d_i, j_i, x_i), i=1,2,\cdots,n$
$t_i$表示時間,
$$d_i=\left\{
\begin{aligned}
& 1 &, \text{if dead} \\
& 0 &, \text{if censored}.
\end{aligned}
\right.
$$
$j_i\in \{1,2,\cdots,M\}$，但對於$d_i=0$即刪失的情況，$j_i$無定義。
4.1 似然函數The likelihood function
\[L=\prod_{i=1}^{n} \lambda_{j_i}(t_i,x_i)^{d_i} S(t_i,x_i):\]
當個體在$t_i$時刻刪失時，體現的是直到$t_i$仍存活的概率:
\[S(t_i,x_i).\]
當個體在$t_i$時刻因原因$j_i$而死亡時，體現的是在$t_i$時刻因原因$j_i$死亡的概率密度$f_{j_i}(t_i,x_i)$:
通過$\lambda_{j_i}(t_i,x_i)$的定義，進行簡單推導可以展開為
\[\lambda_{j_i}(t_i,x_i)\cdot S(t_i,x_i).\]
引入$d_i$，則將這兩個式子統一成同一個形式，得到上面的似然函數
\[L=\prod_{i=1}^{n} \lambda_{j_i}(t_i,x_i)^{d_i} S(t_i,x_i)\]
又由於$S(t_i,x_i)=\prod_j S_j(t_i,x_i)$，則
\begin{eqnarray*}
L & = & \prod_{i=1}^{n} \lambda_{j_i}(t_i,x_i)^{d_i} \prod_j S_j(t_i,x_i) \\
& = & \prod_{i=1}^{n} \lambda_{j_i}(t_i,x_i)^{d_i} \prod_j e^{-\Lambda_j(xt_i,x_i)}
\end{eqnarray*}
用$d_{ij}$表示個體$i$是否因原因$j$而死亡。顯然$d_i=\sum_{j} d_{ij}$.(前面有做設定，個體只能因一種原因死亡)。則有
\[L=\prod_{i=1}^n \prod_{j=1}^M \lambda_j(t_i,x_i)^{d_{ij}} e^{-\Lambda_j(t_i,x_i)}.\]
乘子的順序並不重要，由此還可以引出兩個結論:
1) 總體似然函數是不同failure的似然函數的乘積。這意味着，當不同的$\lambda_j(t,x)$不依賴於相同的參數時，可以對似然逐個求最大值來做估計。
2) 涉及某一類failure的似然函數，和將其他所有failure視為刪失情況下得到的似然函數是完全一樣的。
4.2 Weibull Regression - 韋布爾回歸
假設第j風險函數服從帶Weibull基准風險的比例模型(proportional model with Weibull baseline):
\[\lambda_j(t,x) = \lambda_{j0} e^{x^\mathrm{T} \beta_j}\]
Weibull基准風險為
\[\lambda_{j0}(t) = \lambda_j p_j(\lambda_j t)^{p_j-1}\]
$*$用和前面分析相同的觀點來看，可以同樣把非第j原因造成的死亡視為刪失，對參數$(\lambda_j, p_j,\beta_j)$進行估計。
這里的參數都依賴於j有不同的值，如果有需要的話，對於每種failure，采用不同的x也是可以的。
而如果想讓這些Weibull都有相同的index，比如說$p_j = p_{\forall j}$。那總體的似然函數里會有一些項不能消掉，估計參數也不能用(*)所述的簡化版本了；同樣的，如果讓$\beta$們也共用，也會有這個問題。
(為什么會這樣，文章里沒明講，我覺得是因為如果假設不同的failure里采用的是同一套參數，這就違背了"不同failure之間是獨立的"這個假設。)
4.3 Cox Regression - Cox回歸
如果對$\lambda_{j0}(t)$沒有假設，也可以擬合PH模型(比例模型)。標准的Cox論述導出的偏似然函數如下:
\[L = \prod_{j=1}^m \prod_{i=1}^{k_j} \frac{e^{x^\mathrm{T} _{ji(j)} \beta_j}}{\sum_{k \in R(t_ji)} e^{x^\mathrm{T} _{jk} \beta_j}}\]
其中，$k_j$表示由於原因j導致死亡的次數(相同時間多次的記一次:distinct)，$t_{ji}$表示第i個這樣的時刻，$R(t_{ji})$表示$t_{ji}$時刻的風險集，$i(j)$表示在$t_{ji}$死亡的樣本下標/編號(index)。
如果想將m個基准風險函數都限制為與一個super-baseline成比例，像$$\lambda_{j0}(t) = \lambda_0(t) e^{\gamma_j}.$$
則會得到不同的偏似然(可見K-P中的公式8.15-8.16，不過不知道具體指的是什么，文章沒給出參考文獻或相關鏈接)
4.4 Piece-wise Exponential Survival - 分段指數生存模型
根據Holford，或是Laird和Olivier的論述，用拐點(breakpoints) $0 = \tau_1<\tau_2<\cdots<\tau_k = \infty$定義各個區間，假設j類failure的基准風險服從階梯函數，在各區間中為常數值:
\[\lambda_{j0} = \lambda_{jk}, \mbox{for } t \in [\tau_k, \tau_{k+1}).\]
則似然函數中對應j類failure的因子，和，假設在區間$k$因$j$類failure死亡的人數服從以下式為均值的泊松分布而得到的泊松似然的核，是一樣的:
\[\mu_{ijk} = E_{ik} \lambda_{jk} e^{x_i^\mathrm{T} \beta_j}\]
$E_{ik}$指的是時間處在區間$k$中時，帶有協變量$x_i$的人群里暴露在風險中的數量(total exposure of people with covariates $x_i$ in interval k)
(個體處於風險時，是同時含有因為任何一個原因死亡的風險的，這里$E_{ik}$也即暴露於風險的數量，是不區分原因的(is not cause-specific))
因此，我們可以用 以由於各個原因造成的死亡數作為結果變量、以暴露於所有原因的(數量)為偏置的一系列泊松回歸來擬合競爭風險模型。
(Thus, we can fit competing risk models by running a series of Poisson regressions where we treat the number of deaths due to each cause as the outcome and the exposure to all causes as the offset.)
這個模型的一個優點在於考慮到了：在給定個體因某些原因在時間t死亡的條件下，因原因j而死亡的條件概率。
在給定存活到恰好在t時之前，因原因j而死亡的概率為:
\[\lambda_j(t,x) = \lambda_{j0}(t)e^{x^\mathrm{T} \beta_j} = e^{\alpha_{jk} + x^\mathrm{T} \beta_j},\]
也就是說，$\alpha_{jk} = \log \lambda_{jk}$是在區間k時原因j對應的對數基准風險。
此時在該瞬間的總體風險為
\[\lambda(t,x) = \sum_{j=1}^m \lambda_j(t,x) = \sum_{j=1}^m e^{\alpha_{jk}+x^\mathrm{T} \beta_j}.\]
然后可以得到上述的條件概率:
\[\pi_{jk} = \frac{e^{\alpha_{jk}+x^\mathrm{T}\beta_j}}{\sum_{r=1}^m e^{\alpha_{rk}+x^\mathrm{T} \beta_r}},\]
可以看到這是服從一個以相同的$\beta_j$作為參數的多項Logit模型來作為競爭風險模型!
這意味着我們可以把競爭模型的分析分為兩部分，用
1. 一個用於得到整體風險的標准風險模型
2. 一個關於死亡原因的多項Logit模型
這樣得到的結果，會和直接對每類failure分別擬合泊松模型得到的結果完全一樣。
5. 可識別問題 The Identification Problem
用$T_1, T_2, \cdots, T_m$ 表示個體因$1,2,\cdots,m$原因而死亡的時間。
但顯然對於每一個個體，我們都只能觀察到一個死亡時間點，最早的那一個。
嚴格來說，$T$和$J$都是隨機變量:
\[T = min\{T_1, T_2, \cdots, T_m\}\]
\[J = \{j:T_j\leq T_k \forall k\}\]
5.1 - 5.3 多元生存函數及邊緣分布函數 Multivariate and Maginal Survival
引入聯合生存函數，也稱為Multiple decrement function
\[S_M(t_1, t_2, \cdots, t_m) = Pr\{T_1\geq t_1, T_2 \geq t_2, \cdots, T_m \geq t_m\}.\]
則有
\[S(t) = S_M(t, t, \cdots, t).\]
我們看到，$S(t)$意外地被定義得很妥當(incidentially, S(t) is well defined.)
而case-specific hazards可以從$S_M$的對數偏導入手定義:
\begin{eqnarray*}
\lambda_j(t) & = & \lim_{\Delta t \to 0} \frac{Pr\{t \leq T_j < t+ \Delta t|T>t\}}{\Delta t} \\
& = & - \frac{\partial}{\partial t} \log S_M(t, t, \cdots, t).
\end{eqnarray*}
由於\[Pr\{t\leq T_j < t+\Delta t|T>t\} = Pr\{t\leq t+\Delta t, J=j|T>t\},\]
所以$\lambda_j(t)$與前面部分得到的結果是一樣的。
*要注意的是，由數據對應得到的似然函數只依賴於$\lambda_j(t)$，則只有它們或其函數才能被估計。其他的都是不可識別的(無法估計出來)。
所以，隱藏的各個failure的時間，其邊緣分布是不可識別的:
\[S_j^\ast (t) = Pr\{T_j\geq t\} = S_M(0,\cdots,0, t,0,\cdots,0).\]
$t$出現在第$j$項。
依次可以定義邊緣風險函數:
\[\lambda^\ast (t) = \frac{\partial}{\partial t} \log S_j^\ast (t) \]
一般來說，\[\lambda_j^\ast (t) \neq \lambda_j(t)\]
但如果$T_j$之間互相獨立，則有
\[S_M(t_1,\cdots, t_m) = \prod_{j=1}^m S_j^\ast (t_j).\]
且\[S_j^\ast (t) = S_j(t)\]
此時會有\[\lambda_j^\ast (t) = \lambda_j(t)\]
但是，$S_M$無法被估計，所以這里的獨立性也是沒有辦法去檢驗的。
文章最后舉了一個例子，可能有人會覺得可以通過$\lambda_j(t)$的形式來判斷是否各分量之間是獨立的，但事實上，在一個沒有獨立性假設狀況下推導出的$\lambda_j(t)$，在假設獨立性存在時，仍然可以得到一個不同的$S(t,\cdots,t)$ ，而在這個$S(t,\cdots,t)$中，各分量間確實是獨立的。這意味着獨立性存不存在的兩種情況，可能是會推導出同一個$\lambda_j(t)$的表達式的，試圖由此判斷是不對的。
具體的例子如下:
如果能同時觀察到一個以上的"死亡"時間，那這個識別問題也許是可以解決的。但從定義上看可以知道，一般來說是不可行的。在Panel analysis中，把死亡和消耗都認為是competind risk的一種，則是可以的，因為這兩者可以同時被觀察到。(這個例外沒有接觸，不是很懂具體的情景是什么樣的)
文章最后還有一點討論，關於獨立性(大致結論是 很難 )以及凈概率的定義(計算某類failure的概率時是否考慮了其他failure)。