参考资料:

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)。