生存分析（2）

接上篇文章生存分析（1），本文将进一步介绍生存分析中有关生存率的具体算法及其应用

Kaplan-Meier法（K-M法、乘积极限法）

Kaplan-Meier法由Kaplan和Meier于1958年提出，直接用概率乘法定理估计生存率，故称乘积极限法（product-limit method），是一种非参数法。
这里写图片描述

1.计算方法

a. 将样本生存时间T从小到大排列成如表第1栏。若遇到非删失值和删失值相同时，非截删失排在前面。
b. 列出与T相应的死亡人数d，如表第2栏
c. 列出期初病例数n,如表第3栏，即生存期为某时间时尚存活的病例数
d. 计算活过各时点的生存率P（T>t），计算公式为

P (T > t) = \prod p = \prod n - d n

$P(T>t)=\prod p=\prod {\frac{n-d}{n}}$

2.标准误计算

均数的标准误
为了表示均数的抽样误差大小如何，用的一种指标称为均数的标准误。我们以样本均数为变量，求出它们的标准差即可表示其变异程度，所以将样本均数这“标准差”定名为均数的标准误，简称标准误，以区别于通常所说的标准差。标准差表示个体值的散布情形，而标准误则说明样本均数的参差情况，两者不能混淆。
（具体可参考这里的介绍）

对于K-M法，标准误的计算方法有两种：
<1> $Sp(T>t)=P(T>t)\sqrt{ \sum{\frac{d}{(n-d)n}}}$

<1> $Sp(T>t)=P(T>t)\sqrt{\frac{1-P(t>t)}{n-d}}$

例数较多时，两法计算结果相同，但例数逐渐减少，法1的结果偏小，法二的结果偏大。

3.生存率的可信区间

利用正态近似原理，估计总体生存率的可信区间,如95%置信度

可 信 区 间 = P (T > t) \pm 1.96 S p (T > t)

$可信区间=P(T>t)\pm 1.96Sp(T>t)$
上表中存活时间大于30天，生存率的95%可信区间为

P (T > 30) \pm 1.96 \times 0.1107 = (0.4646, 0.8486)

$P(T>30)\pm 1.96 \times 0.1107=(0.4646,0.8486)$

4.单因素分析（log-rank test）

对数秩检验（log-rank test）
属于非参数检验，用于比较两组或多组生存曲线或生存时间是否相同
检验统计量为卡方 $\chi^2$
自由度=组数-1

χ 2 = \sum ( A - T ) 2 T

$\chi^2=\sum \frac{(A-T)^2}{T}$
其中A为观察死亡数，T为理论死亡数。当有T<5时，用下式进行校正

χ 2 = \sum ( | A - T | - 0.5 ) 2 T

$\chi^2=\sum \frac{(|A-T|-0.5)^2}{T}$
计算完卡方值，查表得到P值，可得到推断结论
具体而言，首先将数据按如下形式组织：

这里写图片描述

将A、B两组的生存天数混在一起从小到大排序放在第二列，第一列是对应的组别，其他按表中给的填入
如此便可以得到A、B两组各自的合计理论死亡数，和实际死亡数（注意删失数据不参与计算），分别带入公式计算卡方即可

K-M方法提供三种假设检验分别是
Log-rank检验、Breslow检验和Tarone Ware检验
三者都是构造卡方检验量，具体比较如下：
这里写图片描述
更为具体的内容可参考这儿

当数据量较小（n<=50）且不含删失数据时，可以选择Wilcoxon 秩和检验，检验效果更好。

Wilcoxon 秩和检验
核心思想：如果两个样本来自相同的整体，那么秩将大约均匀的分布在两个样本中（秩：将样本从小到大排序，排名即为秩）。否则，则有一个样本获得较小的秩和，另一个获得较大的秩和。

计算：
设两个独立样本为：第一个样本x的样本容量为n1，第二个样本y的样本容量为n2，在容量为n1+n2的混合样本（样本x、y之和）中，x样本的秩和为 $W_x$ ，y样本的秩和为 $W_y$ ，且有

W x + W y = 1 + 2 + \dots + n = n ( n + 1 ) 2

$W_x+W_y=1+2+\cdots+n=\frac{n(n+1)}{2}$
我们定义

W1=Wx−n1(n1+1)2=x统一秩和−x原秩和 $W_1=W_x-\frac{n_1(n_1+1)}{2}=x统一秩和-x原秩和$

W2=Wy−n2(n2+1)2=y统一秩和−y原秩和 $W_2=W_y-\frac{n_2(n_2+1)}{2}=y统一秩和-y原秩和$

可以知道样本x在混合之后的秩和最小也是原秩和，即
$min(W_x)=\frac{n_1(n_1+1)}{2}$
对样本y同理，即
$min(W_y)=\frac{n_2(n_2+1)}{2}$

根据 $W_x$ 和 $W_y$ 的关系可知，最大秩和为
$max(W_x)=\frac{n(n+1)}{2}-\frac{n_2(n_2+1)}{2}$
$max(W_y)=\frac{n(n+1)}{2}-\frac{n_1(n_1+1)}{2}$

因此 $W_1$ 和 $W_2$ 的取值范围均为：
$[0,\frac{n(n+1)}{2}-\frac{n_1(n_1+1)}{2}-\frac{n_2(n_2+1)}{2}]=[0,n_1n_2]$

接下来我们进行假设检验。

假设： x，y样本来自相同总体
当原假设为真时，所有的 $x_i$ 和 $y_i$ 相当于从同一总体中抽得的独立随机样本， $x_i$ 和 $y_i$ 构成可分辨的排列情况，可看成一排n个球随机地指定 $n_1$ 个为x球另 $n_2$ 个为y球，共有 $C_n^{n_1}$ 种可能，而且它们是等可能的。基于这样分析，在原假设为真的条件下不难求出 $W_1$ 和 $W_2$ 的概率分布，显然它们的分布还是相同的，这个分布称为样本大小为 $n_1$ 和 $n_2$ 的Mann-Whitney-Wilcoxon分布

一个比较实际的方法是，对于每个样本数大于等于8的大样本来说，我们可以采用标准正态分布Z来近似检验。
因为 $W_1$ 的中心点为 $\frac {n_1n_2}{2}$ ，所以 $W_x$ 的中心点为

μ = n 1 n 2 2 + n 1 ( n 1 + 1 ) 2 = n 1 ( n 1 + n 2 + 1 ) 2

$\mu=\frac{n_1n_2}{2}+\frac{n_1(n_1+1)}{2}=\frac{n_1(n_1+n_2+1)}{2}$
W_x的方差

σ2 $\sigma^2$ 从数学上可推导出

σ 2 = n 1 n 2 ( n 1 + n 2 + 1 ) 12

$\sigma^2=\frac{n_1n_2(n_1+n_2+1)}{12}$

如果样本中存在结，将影响公式中的方差（结：即相同的数据，此时秩会被平分）
按结值调整方差的公式为：

σ 2 = n 1 n 2 ( n 1 + n 2 + 1 ) 12 - n 1 n 2 ( \sum ( τ 3 j - τ j ) 12 ( n 1 + n 2 ) ( n 1 + n 2 - 1 )

$\sigma^2=\frac{n_1n_2(n_1+n_2+1)}{12}-\frac{n_1n_2(\sum(\tau_j^3-\tau_j)}{12(n_1+n_2)(n_1+n_2-1)}$
其中

τj $\tau_j$ 为第

j $j$ 个结的个数。结值的存在将使原方差变小，这是一个显然正确的事实。标准化后

Wx $W_x$ 为

z = W x - μ \pm 0.5 σ = W x - n 1 ( n 1 + n 2 + 1 ) 2 \pm 0.5 n 1 n 2 ( n 1 + n 2 + 1 ) 12 - n 1 n 2 \sum ( τ 3 j - τ ) 12 ( n 1 + n 2 ) ( n 1 + n 2 - 1 ) - - - - - - - - - - - - - - - - - - - - - - - - - \sqrt \sim N (0, 1)

$z=\frac{W_x-\mu\pm0.5}{\sigma}=\frac{W_x-\frac{n_1(n_1+n_2+1)}{2}\pm0.5}{\sqrt{\frac{n_1n_2(n_1+n_2+1)}{12}-\frac{n_1n_2\sum(\tau_j^3-\tau)}{12(n_1+n_2)(n_1+n_2-1)}}}\sim N(0,1)$
其中0.5是为了对离散变量进行连续性修正，对于

Wx−μ>0 $W_x-\mu>0$ 减0.5修正，反之加0.5修正。

算例
x组：11 15 10 18 11 20 24 22 25
y组：13 14 10 8 16 9 17 21
将二者统一如下：

这里写图片描述