生存分析的统计知识小总结

名词解释

起始事件：生存分析中，反应随访开始特征的事件，如疾病确诊、某种治疗开始等，起点一般定义为调查的开始
终点事件：生存分析中，反应生存研究终点的事件，如患者死亡、疾病复发和疾病转移等，这个终点是由研究者定义的
删失事件：发生了非终点事件导致研究终止的事件都叫删失事件，删失有三类，如果目标事件在目前知晓时间之前发生，则为左删失，如果目标事件在某一区间发生，成为区间删失，如果目标事件在目前知晓时间之后发生，则为右删失，下面用直白的语言总结一下常见的删失事件：
- 左删失事件：病人死于最后一次随访之前，但不清楚具体的死亡时间
- 右删失事件：研究截止日期时，病人还活着、病人在最后一次记录后失去联系、病人中途退出研究等等
生存时间：从起点时间到终点事件之间经历的时间跨度为生存时间
中位生存时间：使得刚好有一半的个体的存活时间大于该时间的时间点
竞争风险：当研究对象会发生多结局事件，并且一个结局事件的发生会阻碍其他结局事件发生时（或者这些结局事件间互斥），便称这些结局事件间存在竞争风险。在生存分析中，只有以全死因死亡为结局时，才不存在竞争风险。以死因别死亡（如心血管疾病死亡）为结局或者以疾病发病为结局，均存在竞争风险。

常用公式

生存概率： $p = \frac {活过该时间点人数} {初期人数}$
累积生存率： $S(t, \boldsymbol{X}) = P(T>t, \boldsymbol{X})$ ，表示具有协变量 $\boldsymbol{X} = (X_1, X_2, \dots, X_m)'$ 的观察对象的生存时间 $T$ 大于某时刻 $t$ 的概率， $X_i$ 表示一种因素，也称为生存函数
死亡概率： $q = 1 - p$
累积死亡率： $F(t, \boldsymbol{X}) = P(T \le t, \boldsymbol{X})$ ，表示具有协变量 $\boldsymbol{X} = (X_1, X_2, \dots, X_m)'$ 的观察对象的生存时间 $T$ 不大于某时刻 $t$ 的概率，也称为死亡函数，数值上等于 $1-S$
风险函数： $h(t, \boldsymbol{X}) = \lim_{\Delta t \to 0} \frac {P(t<T\le\Delta t | T \ge t, \boldsymbol{X})} {\Delta t} = \frac {f(t, \boldsymbol{X})} {S(t, \boldsymbol{X})}$ ，表示具有协变量 $\boldsymbol{X} = (X_1, X_2, \dots, X_m)'$ ，且生存时间已达到 $t$ 的观察对象在时刻 $t$ 的瞬时死亡率，实际上是一个条件概率
累积风险函数： $H(t, \boldsymbol{X}) = -lnS(t, \boldsymbol{X})$ 或 $S(t, \boldsymbol{X}) = exp(-H(t, \boldsymbol{X}))$

分析方法

Kaplan-Meier法#

　　Kaplan-Meier法由Kaplan和Meier在1958年提出，该方法也叫乘积极限法，是一种估计生存率及其标准误，计算可信区间的方法，根据计算结果可以绘制表示生存率变化的图，叫做KM曲线，是目前展示生存分析结果的最常用的图表类型，以下的内容逐步展示KM算法，假设我们有两组数据，生存时间分别是：

甲组 (n = 23)：1 3 5 5 5 6 6 6 7 8 10 10 14+ 17 19+ 20+ 22+ 26+ 31+ 34 34+ 44 59
乙组 (n = 20)：1 1 2 3 3 4 4 4 6 6 8 9 9 10 11 12 13 15 17 18

排序与合并#

　　拿到患者的生存时间数据后，先把所有患者的生存时间进行从小到大排序，并列出所有的生存时间以及对应的人数，如果两个或多个患者具有相同的生存时间但性质不同，则算作两个时间，这一过程主要统计以下表格（以甲组为例），数字表示生存时间，数字的个数表示死亡数，如果有删失，就记为删失：

序号 $i$	时间 $t_i$	死亡 $d_i$	删失 $c_i$	期初 $n_i$
1	1	1	0	23
2	3	1	0	22
3	5	3	0	21
4	6	3	0	18
5	7	1	0	15
6	8	1	0	14
7	10	2	0	13
8	14	0	1	11
9	17	1	0	10
10	19	0	1	9
11	20	0	1	8
12	22	0	1	7
13	26	0	1	6
14	31	0	1	5
15	34	1	0	4
16	34	0	1	3
17	44	1	0	2
18	59	1	0	1

生存概率、死亡概率、累积生存率和标准误#

生存概率： $p_i = (n_i-q_i-c_i) / (n_i)$
死亡概率： $q_i = 1 - p_i$
累积生存率： $\hat{S} (t_i) = \prod_{n=1}^n p_i$
标准误： $SE(\hat{S} (t_i)) = \hat{S}(t_i) \sqrt{\sum_{j=1}^i \frac {d_j} {n_j(n_j - d_j)}}$

　　这一过程后，表格如下：

序号 $i$	时间 $t_i$	死亡 $d_i$	删失 $c_i$	期初 $n_i$	生存概率 $p_i$	死亡概率 $q_i$	生存率 $\hat{S}(t_i)$	标准误 $SE[\hat{S}(t_i)]$
1	1	1	0	23	0.956521739	0.043478261	0.95652	0.042523
2	3	1	0	22	0.954545455	0.045454545	0.91304	0.058753
3	5	3	0	21	0.857142857	0.142857143	0.78261	0.086006
4	6	3	0	18	0.833333333	0.166666667	0.65217	0.099311
5	7	1	0	15	0.933333333	0.066666667	0.6087	0.101764
6	8	1	0	14	0.928571429	0.071428571	0.56522	0.103367
7	10	2	0	13	0.846153846	0.153846154	0.47826	0.104159
8	14	0	1	11	1	0	0.47826	0.104159
9	17	1	0	10	0.9	0.1	0.43043	0.104146
10	19	0	1	9	1	0	0.43043	0.104146
11	20	0	1	8	1	0	0.43043	0.104146
12	22	0	1	7	1	0	0.43043	0.104146
13	26	0	1	6	1	0	0.43043	0.104146
14	31	0	1	5	1	0	0.43043	0.104146
15	34	1	0	4	0.75	0.25	0.32283	0.121597
16	34	0	1	3	1	0	0.32283	0.121597
17	44	1	0	2	0.5	0.5	0.16141	0.129319
18	59	1	0	1	0	1	0	0

置信区间的计算#

\hat{S} (t_{i}) \pm z_{α / 2} \cdot S E [\hat{S} (t_{i})]

$\hat{S} (t_i) \pm z_{\alpha/2} \cdot SE[\hat{S}(t_i)]$

　　在 $\alpha = 0.05$ 时， $z_{\alpha/2} = 1.96$ ，所以置信区间为生存率加减1.96倍的标准误，但是这种算法有可能使置信区间超出 $[0,1]$ 的范围，不合常理，因此我们可以使用Kalbfleisch和Prentice提出的新算法来计算对数对数生存率的可信区间，这个算法可以表示为：

S E [l n [- l n \hat{S} (t_{i})]] = \frac{S E [\hat{S} (t_{i})]}{| l n \hat{S} (t_{i}) | \cdot \hat{S} (t_{i})}

$SE[ln[-ln\hat{S}(t_i)]] = \frac {SE[\hat{S}(t_i)]} {|ln\hat{S}(t_i)| \cdot \hat{S}(t_i)}$

　　对应的95%置信区间为：

l n [- l n \hat{S} (t_{i})] \pm 1.96 * \frac{S E [\hat{S} (t_{i})]}{| l n \hat{S} (t_{i}) | \cdot \hat{S} (t_{i})}

$ln[-ln\hat{S}(t_i)] \pm 1.96 * \frac {SE[\hat{S}(t_i)]} {|ln\hat{S}(t_i)| \cdot \hat{S}(t_i)}$

　　转换为生存率的置信区间就是：

e x p {- e x p [l n [- l n \hat{S} (t_{i})] \pm 1.96 * \frac{S E [\hat{S} (t_{i})]}{| l n \hat{S} (t_{i}) | \cdot \hat{S} (t_{i})}]}

$exp\{ -exp[ ln[ -ln \hat{S} (t_i) ] \pm 1.96 * \frac {SE[\hat{S}(t_i)]} {|ln\hat{S}(t_i)| \cdot \hat{S}(t_i)} ] \}$

组间生存时间差异的检验#

Log-rank检验#

　　Log-rank检验是目前最常用的检验方法，此检验的假设是：

$H_0$ ： $S_1(t) = S_2(t)$
$H_1$ ： $S_1(t) \ne S_2(t)$

　　对应的步骤如下：

按照上文中的方式将数据进行排序和整理，需要以下数据：
- 序号： $i$
- 生存时间： $t_i$
- 各组的死亡数： $d_{ki}$
- 各组的删失数： $c_{ki}$
- 各组的期初样本： $n_{ki}$
- 各组的理论死亡数： $T_{ki}$ ，这个数字是用所有组的死亡数之和除以所有组的期初样本之和得到死亡比率，然后分别乘以每组的期初样本数： $T_{ki} = n_{ki}\frac {\sum_k d_{ki}} {\sum_k n_{ki}}$
计算方差估计值： $V_{ki} = \frac {n_{ki}} {n_i}(1 - \frac {n_{ki}} {n_i})(\frac {n_i - d} {n_i - 1})d_i$
计算 $\chi^2$ 值： $\chi^2 = \frac {(\sum d_{ki} - \sum T_{ki})^2} {\sum V_{ki}}$ ，自由度为组数减一，使用R中的 1 - pchisq(n, df) 可计算 $p$ 值

Breslow检验#

　　Breslow检验与Log-rank检验的不同之处是Breslow检验在计算卡方统计量的时候加入了权重：

χ^{2} = \frac{(\sum w_{i} d_{k i} - \sum w_{i} T_{k i})^{2}}{\sum w_{i}^{2} V_{k i}}

$\chi^2 = \frac {(\sum w_i d_{ki} - \sum w_i T_{ki})^2} {\sum w^2_i V_{ki}}$

　　其中，权重 $w_i = \sum_k n_{ki}$ 即时间点 i 各组期初样本数的总和，考虑到随着 i 增大，期初样本会减少，因此Breslow检验对前期的差别更加明显，相比之下，Log-rank检验对于全程的权重是相同的

绘制生存曲线#

　　生存曲线指的是KM曲线，以生存时间为横轴，生存率为纵轴

计算中位生存时间#

　　绘制出生存曲线后，可以先在纵轴的50%处绘制一条垂直于纵轴的直线，在直线与生存曲线的交点绘制垂直于横轴的直线，该直线与横轴的交点就是中位生存时间

　　如果要获得具体的数值，可以使用线性内插法计算中位生存时间，首先找出两个生存率 $S(t_{i-1})$ 和 $S(t_i)$ ，使得 $S(t_{i-1}) > 0.5$ 且 $S(t_{i}) < 0.5$ ，然后列一个方程求出中位生存时间 $x$ ：

\frac{t_{i - 1} - t_{i}}{t_{i - 1} - x} = \frac{S (t_{i - 1}) - S (t_{i})}{S (t_{i - 1}) - 0.50}

$\frac {t_{i-1} - t_i} {t_{i-1} - x} = \frac {S(t_{i-1}) - S(t_{i})} {S(t_{i-1}) - 0.50}$

Cox比例风险回归#

　　生存分析的主要目的在于研究变量与累积生存率之间的关系，当累积生存率受很多因素影响时，传统的方法是使用多元回归检查每一个因素对因变量的影响。但由于生存分析研究中的数据包含删失数据。且时间变量 $t$ 通常不满足正态分布和方差齐性的要求，这就使得用一般的回归方法研究上述关系比较困难，因此D.R.Cox提出了Cox比例风险回归模型，

Cox模型#

　　与传统生存分析的区别在于选择风险函数作为因变量，该模型可以描述为：

\begin{matrix} (1) & h (t, X) = h_{0} (t) e x p (β^{'} X) = h_{0} (t) e x p (β_{1} X_{1} + β_{2} X_{2} + \dots + β_{m} X_{m}) \end{matrix}

$\begin{equation} h(t, \boldsymbol{X}) = h_0(t)exp(\boldsymbol{\beta'X}) = h_0(t)exp(\beta_1X_1 + \beta_2X_2 + \dots + \beta_mX_m) \end{equation}$

　　其中， $h(t, \boldsymbol{X})$ 是具有协变量 $\boldsymbol{X}$ 的个体在时刻 $t$ 时的风险函数； $t$ 是生存时间； $\boldsymbol{X} = (X_1, X_2, \dots, X_m)'$ 是可能影响生存时间的有关因素，或称协变量，这些变量可以是定量的，也可以是定性的，并且在整个观察期间不随时间而变化； $h_0(t)$ 表示当所有协变量为0时的风险函数，称为基线风险函数（baseline hazard function）； $\boldsymbol{\beta} = (\beta_1, \beta_2, \dots, \beta_m)'$ 是Cox模型的回归系数，是一组待估计的回归参数。

　　用生存率表示的Cox模型为：

\begin{matrix} (2) & S (t, X) = S_{0} (t)^{e x p β^{'} X} = S_{0} (t)^{e x p (β_{1} X_{1} + β_{2} X_{2} + \dots + β_{m} X_{m})} \end{matrix}

$\begin{equation} S(t, \boldsymbol{X}) = S_0(t)^{exp{\boldsymbol{\beta'X}}} = S_0(t)^{exp(\beta_1X_1+ \beta_2X_2 + \dots + \beta_mX_m)} \end{equation}$

　　上式中， $S(t, \boldsymbol{X})$ 是具有协变量 $\boldsymbol{X}$ 的个体在时刻 $t$ 的生存率， $S_0(t)$ 是在 $t$ 时刻的基线生存率，其他符号不变

　　由于Cox回归模型对 $h_0(t)$ 未做任何假定，不需要服从特定的分布，具有非参数的特点，而指数部分 $exp(\boldsymbol{\beta'X})$ 具有参数模型的形式，因此整个模型是半参数模型

Cox回归的假设#

　　使用Cox回归模型之前，需要判断所研究的数据是否符合以下假设：

比例风险假定：各危险因素的效应不随时间而变化
对数线性假定：模型中的协变量与对数风险比应该呈线性关系

回归系数的估计#

　　使用偏似然函数和最大似然估计可以获得Cox模型中各协变量的回归系数，偏似然函数的计算公式为：

\begin{matrix} (3) & L = q_{1} q_{2} \dots q_{i} \dots q_{k} = \prod_{i = 1}^{k} q_{i} = \prod_{i = 1}^{k} \frac{h (t_{i})}{\sum_{j = i}^{n} h_{j} (t)} = \prod_{i = 1}^{k} \frac{e x p (β_{1} X_{i 1} + β_{2} X_{i 2} + \dots + β_{m} X_{i m})}{\sum_{S \in R (t_{i})} e x p (β_{1} X_{s 1} + β_{2} X_{s 2} + \dots + β_{m} X_{s m})} \end{matrix}

$\begin{equation} L = q_1q_2 \dots q_i \dots q_k = \prod_{i = 1}^{k} q_i = \prod_{i = 1}^{k} \frac {h(t_i)} {\sum_{j=i}^{n} h_j(t)} = \prod_{i = 1}^{k} \frac {exp(\beta_1X_{i1} + \beta_2X_{i2} + \dots + \beta_mX_{im})} {\sum_{S \in R(t_i)} exp(\beta_1X_{s1} + \beta_2X_{s2} + \dots + \beta_mX_{sm})} \end{equation}$

　　其中， $q_i$ 是第 $i$ 个死亡时间点的条件死亡概率，它的分子部分是第 $i$ 个个体在 $t_i$ 死亡时间点的风险函数 $h(t_i)$ ，它的分母部分是处于风险的个体（生存时间 $T \ge t_i$ 的所有个体）的风险函数之和，对函数 $L$ 取对数，得到对数似然函数 $lnL$ ，对 $lnL$ 求关于 $\beta_j \space (j=1,2, \dots, m)$ 的一阶偏导数，并对 $\frac {\partial{lnL}} {\partial{\beta_j}} = 0$ 求解，便可获得 $\beta_j$ 的最大似然函数估计值 $b_j$

回归系数的检验#

似然比检验
Wald检验
计分检验

相对危险度和风险比#

　　当我们想判断众多协变量中每一个协变量分别对生存的影响，我们就需要计算相对危险度 (RR) 或风险比 (HR)，表示观测值每增加一个单位时引起的累积生存率的变化，将公式（1）进行变形，可得：

\begin{matrix} (4) & l n [\frac{h (t, X)}{h_{0} (t)}] = l n R R = β_{1} X_{1} + β_{2} X_{2} + \dots + β_{m} X_{m} \end{matrix}

$\begin{equation} ln[\frac {h(t, \boldsymbol{X})} {h_0(t)}] = lnRR = \beta_1X_1 + \beta_2X_2 + \dots + \beta_mX_m \end{equation}$

　　假设有两个分组，对照组和观察组， $X_j$ 表示对照组中的各因素的取值， $X_i$ 表示观察组中的各因素取值，那么，观察组相对于对照组的相对危险度就可以通过下面的公式求得：

\begin{matrix} (5) & R R = \frac{h (t, X_{i})}{h (t, X_{j})} = e x p [β (X_{i} - X_{j})], i, j = 1, 2, \dots, m \end{matrix}

$\begin{equation} RR = \frac {h(t, \boldsymbol{X_i})} {h(t, \boldsymbol{X_j})} = exp[\boldsymbol{\beta(X_i - X_j})], \space i,j = 1,2,\dots, m \end{equation}$

　　或

\begin{matrix} (6) & R R_{j} = e x p [β_{j} (X_{j} - X_{j}^{*})] \end{matrix}

$\begin{equation} RR_j = exp[\beta_j(X_j - X_j^*)] \end{equation}$

　　所以偏回归系数 $\beta_j$ 的含义就是在其他协变量不变的情况下，协变量 X_j 每增加一个单位时引起的相对危险度的自然对数的该变量，根据以上公式，可以总结协变量大小和危险因素的关系如下：

$\beta_j > 0$ 时，则 $RR_j > 1$ ，则 $X_j$ 为危险因素
$\beta_j = 0$ 时，则 $RR_j = 1$ ，则 $X_j$ 为无关因素
$\beta_j < 0$ 时，则 $RR_j < 1$ ，则 $X_j$ 为保护因素

预后系数#

　　Cox模型的线性部分 $\beta_1X_1 + \beta_2X_2 + \dots + \beta_mX_m$ 与风险函数 $h(t)$ 成正比，因此这个线性部分反映了一个个体的预后，有人将这部分内容定义为预后指数 (prognosis index, PI)，用来表征患者的预后风险，如果对各变量进行标准化转换后再拟合Cox模型，则可得到标准化后的预后指数，此时，标准化预后风险 $PI'$ 与0的比较可以表示患者死亡风险与平均水平的关系