生存分析（2）

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

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

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

1.计算方法

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

2.标准误计算

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

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

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

3.生存率的可信区间

$$可信区间=P(T>t)\pm 1.96Sp(T>t)$$

上表中存活时间大于30天，生存率的95%可信区间为

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

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

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

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

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

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

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

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

$$W_x+W_y=1+2+\cdots+n=\frac{n(n+1)}{2}$$

我们定义
$W_1=W_x-\frac{n_1(n_1+1)}{2}=x统一秩和-x原秩和$
$W_2=W_y-\frac{n_2(n_2+1)}{2}=y统一秩和-y原秩和$

根据$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$的中心点为

$$\mu=\frac{n_1n_2}{2}+\frac{n_1(n_1+1)}{2}=\frac{n_1(n_1+n_2+1)}{2}$$
W_x的方差 $\sigma^2$从数学上可推导出

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

$$\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)}$$
其中 $\tau_j$为第 $j$个结的个数。结值的存在将使原方差变小，这是一个显然正确的事实。标准化后 $W_x$为

$$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是为了对离散变量进行连续性修正，对于 $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
将二者统一如下：

可以看出$n_1=9$，$n_2$=8，$W_x$=96.5，$W_y$=56.5.
H0:两个样本的分布是相同的。标准分布z值的计算结果为：

$$z=\frac{96.5-\frac{9(9+8+1)}{2}-0.5}{\sqrt{\frac{9\times8(9+8+1)}{12}-\frac{9\times8(8-2+8-2)}{12(9+8)(9+9-1)}}}=1.44515$$
如果设定显著水平 $\alpha=0.05$，我们知道标准正态分布在0.05显著水平时，上临界值为1.645，下下临界值为-1.645，由1.445<1.645，所以不能拒绝原假设。
当然你也可以用第二个样本的秩和 $W_y$来计算标准正态分布，此时要注意公式中的 $n_1$和 $n_2$发生对换。

寿命表法（life table，LT）

1.计算方法

2.标准误计算

$$Sp(T>n)=p(T>n)\sqrt {\sum \frac {q}{pN}}$$

3.生存率可信区间

同K-M法

K-M法与寿命表法比较

1.格式与精确度
K-M法使用患者实际寿命作为分布区间，相对更精确
寿命表法采用人为规定时间段作为分布区间，范围扩大精确度有所下降

2.适用范围
K-M法更适合于样本量较少的数据
寿命表法更适合于样本量较大的数据
（不过，考虑到现在计算能力的强大，一般程度的数据量并不会对计算速度有太大影响）

3.关注点不同
K-M法关注每一个时点的生存率，重视对生存率规律的细致把握，可以利用K-M的结果去研究影响生存率变化（如曲线的突变点）的影响因素。
寿命表法则更重视对生存规律的总体把握（如各年生存率的情况）。

至此，有关生存分析的非参数研究方法K-M与寿命表法就介绍完了。接下来的一篇文章，将介绍含参数的研究方法—Cox比例风险回归模型。