参考资料:
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)。