生存分析第二课:Kaplan-Meier estimate of S(t)

Kaplan-Meier estimate of $S(t)$

---

导言

    生存函数$S(t)$是个体至少生存到时刻$t$的概率。
    对$S(t)$的估计可以通过参数视角,也可以通过非参数视角.
    参数视角估计$S(t)$:先假设一个生存时间$T$服从的分布,然后根据样本数据估计出未知参数,最后可估计出$S(t)=1-F(t)$.
    参数视角估计$S(t)$的优点: 模型明确,即使数据量很少也能识别出模型等。
    参数视角估计$S(t)$的缺点: 如果最初假设的那个模型就是错的，那么后续的估计与分析就是徒劳的！
    非参数视角估计$S(t)$: 即使我们没有识别出生存时间$T$的分布,也可以进行估计.
    非参数视角估计$S(t)$的优点: 适用性强,不需要事先假设出明确的模型,出"大错"的概率小。
    非参数视角估计$S(t)$的缺点: 所需样本量较大,样本量小时效果不好,不如参数模型明确.

---

引例: leukæmia(白血病)

    将42名青年分配到实验组和对照组，实验组接受 6-Mercaptopurine (6-巯基嘌呤)的处理,对照组接受 placebo(安慰剂)的处理.
实验结果如下:

6-MP	6, 6, 6, 7, 10, 13, 16, 22, 23, 6+, 9+, 10+, 11+, 17+, 19+,20+, 25+, 32+, 32+, 34+, 35+
Placebo	1, 1, 2, 2, 3, 4, 4, 5, 5, 8, 8, 8, 8, 11, 11, 12, 12, 15, 17, 22, 23
注	数字后面的 + 号表示至少活到这个时刻,只是因为某些原因没有继续被观察到从而删失了

问题1:
    安慰剂组的个体活到以下时刻的比例是多少?

• time 0.0? — 100%
• time 0.9? — 100%
• time 1.0? — 100%
• time 1.1? — 19/21*100≈90.48%
• time 2.0? — ≈90.48%

    总结:如果没有删失情况的存在,那么$\hat{S}$(t)就是一个阶梯函数,每个时间点都对应着一个生存比例。

    但是如果像实验组那样存在删失情形该怎么办呢?

---

Kaplan-Meier’s method

符号引入:
$d(t)$: 时刻$t$ 死亡或失败的个体数;
$q(t)$: 时刻$t$ 右删失的数目;
$n(t^{-})$ :时刻$t$左侧一点点时间时处于风险中的个体数

With right-censored data:

• 公式1— KM estimate of $S(t)$ is $\hat{S}$(t)=$\hat{S}$($t^{-})\hat{p}(T>t|T\ge t)$
• 公式1解读:活过时刻$t$的概率=活过时刻t前一点点的概率×活过时刻t前一点点为条件并活过时刻$t$的概率.
• 其中
    $\hat{p}(T>t|T\ge t)=1$,如果$t$时刻没有任何failure事件发生;
    $\hat{p}(T>t|T\ge t)=\frac{n(t^{-})-d(t)}{n(t^{-})}$,如果$t$时刻没有任何failure事件发生;
• 很明显,KM 估计的$\hat{S}$(t)=$\hat{S}$($t^{-})\hat{p}(T>t|T\ge t)$,对于从时刻$t^{-}$到时刻$t$而言，只有当时刻$t$发生死亡事件,$\hat{p}(T>t|T\ge t)=\frac{n(t^{-})-d(t)}{n(t^{-})}$才会发生从1变为其他数字,$\hat{S}$(t)也才会发生改变.
• 公式2—与KM估计等价的公式:$\hat{S}(t)={\prod }_{t_{(i)\le t}}\frac{n_{(i^{-})}-d_{(i)}}{n_{(i^{-})}}$
• 公式2解读: 活过时刻t的概率=同时满足活过之前每一个时刻=活过之前每一个概率的乘积

---

     按照公式,当时间$t$超过最大观察时间$t_0$时,$\hat{S}(t)$便不再发生形状变化,很有可能还不为0，但这明显与"人都会死""这样的事实相矛盾,这样看来KM 估计似乎不太好，我们该如何处理时间超过最大观察期后的生存函数取值呢?

目前有几种办法处理这种情况:

• Efron (1967) 建议将$t>t_{max}$的$\hat{S}(t)$一律设置为0；
• Gill (1980) 建议将$t>t_{max}$的$\hat{S}(t)$设置为$\hat{S}(t)=\hat{S}(t_{max})$
• Brown et al.(1974)建议将$t>t_{max}$的$\hat{S}(t)$设置为 S ^ ( t ) = e x p { l o g ( S ^ ( t m a x ) ) × t / t m a x } \hat{S}(t)=exp\{log(\hat{S}(t_{max}))×t/t_{max}\} S^(t)=exp{ log(S^(tmax​))×t/tmax​}
• 但事实上,最好的处理方法是不对没有任何信息的事情做任何假设,直接在最大观察点停止绘图.

---

Variance of the KM estimate

• 引言
      对统计学家来说,会自然的想要知道我们的估计有多大的把握-------置信区间
      计算置信区间有一个必须要做的事情就是计算统计量的方差,所以这一部分让我们来计算一下KM估计的方差：

• 推导过程
  根据公式2: V { S ^ ( t ) } = V { ∏ t ( i ) ≤ t n ( i − ) − d ( i ) n ( i − ) } = V { ∏ t ( i ) ≤ t p ^ ( i ) } V\{\hat{S}(t)\}=V\{\prod_{t_{(i)}≤t}{\frac{n_{(i^-)}-d_{(i)}}{n_{(i^-)}}}\}=V\{\prod_{t_{(i)}≤t}{\hat{p}_{(i)}}\} V{ S^(t)}=V{ t(i)​≤t∏​n(i−)​n(i−)​−d(i)​​}=V{ t(i)​≤t∏​p^​(i)​}
      根据统计学常识,计算一系列独立随机变量和的方差是容易的,但是计算连乘的方差是非常困难的,所以我们需要利用取对数变乘为加这个技巧进行转换:
  V { l o g S ^ ( t ) } = V { ∑ t ( i ) ≤ t l o g p ^ ( i ) } = ∑ t ( i ) ≤ t V { l o g p ^ ( i ) } V\{log\hat{S}(t)\}=V\{\sum_{t_{(i)}≤t}{log\hat{p}_{