Censored data：当注射格列宁的小白鼠逃出实验室
对应知识点：CT3-10 Point Estimation中，对截尾和删失数据的极大似然估计方法；
CT4-08 Estimating the lifetime distribution function中，对删失数据的概念和分类的系统介绍。（属于生存分析的内容）通过小白鼠的例子介绍卡普兰·迈耶生存分析法（乘积极限法）。
同时补充了截尾数据和删失数据的对比。
1
我已经有好几天没看到那只小白鼠了。
传言说它逃出实验室是因为要去寻找防治脱发的丹药，也有人说看到它前几天被黑猫警长抓走了，因为躲在图书馆偷看了几本SOA的禁书。世事难料，要知道它原先可是IFoA阵营的忠实粉丝啊。
但准确原因已经无从查考了。我只记得最后一次见到它那天，天气很热。它穿着一件粉红色的T恤，端正地坐在自己的笼子里，时不时地发出“吱吱”的叫声，像极了等着妈妈来接自己放学的小女孩，模样很是乖巧。
“这么热要不给它买杯琪琪家的冰豆浆吧。”我这样想着，就走出了实验室，并给它打开了笼子，“这样它应该也会凉快些。”
其实我认识这只小白鼠的时间并不算长。但对于这只小白鼠，我总有一种熟悉的感觉。
因此我常常给它开个小灶。每次看着它喜庆的吃相，我都会想，也许这就是作为一名饲养员的乐趣吧。
但那天，粗心的我没有注意到，它手心里已经偷偷攥好了一把粉红色的遮阳伞。貌似还有把玩具车的车钥匙。
谁能想到它就这样逃出了实验室呢。
这个不听话的小白鼠。30天的观测期还没结束啊喂！自己还整了个车车给跑了？！
下回再碰到，看我不把它所有的太阳伞都给没收了！
2
故事的起因是新药格列宁的研发。
为测试药效，“精算后花园”实验室对15只小白鼠注射了格列宁。
实验员Jackie对这15只小白鼠进行了为期30天的观测，希望得到小白鼠注射新药后的生存数据（即进行死亡调查）（mortality investigation）。
在这30天里，小白鼠可能死，也可能还活着；死因可能是注射了格列宁，也可能是其他；可能还活着，但逃出了实验室，我们无法再观测到它。
在观测期间内，除了因为注射格列宁而直接引起的死亡以外的所有情况，都会导致出现数据的删失（truncated）或截尾（censored）。
讲到这里我们先把小白鼠的故事放一放，在〔8〕生存分析会继续讲小白鼠们的具体死亡和删失情况。
〔9〕和〔10〕也是关于小白鼠的后续，算是讲了一个完整的故事。如果想先看完故事，可以直接滑到最后。
3
先来讲一下删失。这部分是CT4要求掌握的。
删失（Censoring）
Censored data arises when an individual’s life length is known to occur only in a certain period of time.
也就是说，如果我们想要知道对象存活到的准确的时间点（lifetime），那么事件的失效时间点（failure time）必须要在删失时间的区间以内，也就是要在观察区间以内。
如果失效事件发生的时间不是在观察的时间区间以内，该对象的数据就会以删失数据（censored data）的形式出现。
我们只知道在观测区间以外发生了失效事件，但并不知道具体事件发生的失效时间点的精确值。
和截尾不同的是，数据存在的事实我们是知道的（后面会说截尾的特点，主要是对于观测外的事件的发生并不知情/不考虑）。（不妨思考一下空白数据和0的区别）
以下是删失的几种分类：（注意它们并不是严格互斥的情况）
左删失（Left censoring）
当事件发生时，观测还没有开始，事件发生时间位于观测区间的左侧，那么相应的数据就是左删失的数据。我们只知道事件在观测开始之前就发生，但还是不知道精确的发生时间。
左删失一般发生在：
- 当估计某个确定的年龄的生存函数，但你不知道准确的出生日期时；
- 当估计某个确定的保险期间的生存函数，但你不知道准确的保单起始日期时；
- 当估计某种疾病的持续期间的生存函数，但你不知道准确的开始生病的时间时。
右删失（Right censoring）
在观测区间中，对象并没有发生失效事件，失效事件的发生时间是在我们观测区间的右侧，但我们也不能得到精确的失效时间。
比如在这30天的死亡调查中，当时间超过30天时，也许还是有一部分小白鼠没有死亡，但我们实验的观测和记录已经终止了，那么这些还未死亡的小白鼠的数据就是右删失的数据。我们只知道这些小白鼠的生存时间超过30天。
右删失一般发生在：
- 寿险的投保人退保时；
- 养老金计划的被保险人在退休领取养老金时为健康状态（active）；
- 生死两全保险的保单期满。
区间删失（Interval censoring）
在某一个时间区间内，我们无法知道事件的精确失效时间，只知道失效事件是在这个区间内发生的。
区间删失一般发生在：
- 当你只知道退保的会计年度（calendar year ）时。
- 当估计某个确定的年龄的生存函数，但你只知道死亡发生时的年龄是“x nearest birthday”时。也就是说你只知道死亡发生在【x-0.5，x+0.5】的年龄区间（Rate interval）内；
- 当你知道死亡和出生的会计年度时。因为此时我们有的数据只是“age next birthday at the 1 January prior to death”，我们只知道死亡发生在这个时间区间内。
在精算调查中，右删失是最为常见的。
随机删失（Random censoring）
顾名思义，此时删失是个随机变量，与死亡无关。比如小白鼠因为各种非格列宁引起的原因离开，我们无法继续观测。而这种离开我们事先是不知道的。
随机删失一般发生在：
- 寿险的投保人退保时；
- 人群发生迁移；
- 企业年金计划的成员跳槽到另一家企业。
随机删失是右删失的一种特例。
当删失成为我们考虑的失效事件的第二原因时，就可以建立多重减因模型。这在CT5寿险精算中会详细讲。
Informative and non-informative censoring
如果删失数据不会给出生存时长（余命）任何信息，我们就称这种删失为non-informative censoring.
也就是说，对于non-informative censoring，我们观测的整个群体的死亡率和其中删失的部分的死亡率是一样的。
non-informative censoring的例子：
- 调查期间的结束（因为这对于所有的个体的影响是相同的）。
Informative censoring的分析则比较复杂。（对，你可以把这句话大致理解为Informative censoring只要掌握概念就好）
Informative censoring的例子有：
- 寿险的退保。因为退保的人可能会比没有退保的人健康状态略好，死亡率低一些；
- 养老金计划的被保险人在退休领取养老金时为生病状态（Ill-health）。因为他们的健康状态一般会比留在养老金计划中的人差，死亡率更高。
第一类删失（Type I censoring）
如果删失发生的时刻是事先知道的，我们就称之为第一类删失。
这是随机删失的一种退化的情形。因此，第一类删失也是右删失的一种特例。
第二类删失（Type II censoring）
如果删失发生的数量是事先预定的，我们就称之为第二类删失。
在精算调查中，最常见的是随机删失和第一类删失的组合。
4
CT4对截尾并没有详细讲解，但这里还是补充一下。
截尾/截断（Truncation）
Truncation of survival data occurs when only those individuals whose event time lies within a certain observational window (YL,YR) are observed. An individual whose event time is not in this interval is not observed and no information on this subject is available to the investigator.
也就是说，我们只能从一个固定的时间窗内观测，记录事件在这个时间段内发生的精确时间点。但是对于时间窗以外的时间发生的失效事件，我们是不知情的。只有当事件的失效时间（也就是事件的发生时间，event time）出现在观测区间内，我们才能知道这个事件及其观测数据的存在。
左截尾（Left truncation）
只能观测到一个时间点之后发生的失效事件。左截尾时间点之前发生的失效事件不知情/不关心。
右截尾（Right truncation）
只能观测到一个时间点之前发生的失效事件。右截尾时间点之后发生的失效事件不知情/不关心。
区间截尾（Interval truncation）
在这个固定区间内发生的失效事件不知情，也就是在这个区间内发生的失效事件不能被观测到，或者说不关心或不知情。
截尾（truncated）是所有样本的综合特性，指的是你观察的总体是有偏的。比如你的寿命样本全部来自于退休中心，但是退休中心只有60岁以上才能加入，所以你的所有样本都是t≥60，t＜60的数据被排除掉了，是左截尾（注意此处与删失不同）；要是你的所有样本只是80岁以下死亡的老人的寿命，那就面临右截尾。
5
删失和截尾最主要的区别是什么呢？
观测者对于删失数据的发生是知情的，只是仍然不知道对应这些删失数据，其失效事件发生的精确时间。
而观测者对于截尾数据的发生是不知情或不关心的。
因此构造似然函数时，截尾数据要引入条件分布，而删失数据要考虑累积概率分布。
6
既然讲到了构造似然函数，我们就来看看对于删失数据和截尾数据，我们分别应该如何进行极大似然估计。这是CT3-10 Point Estimation的内容。
极大似然估计法的基本步骤
我们知道，极大似然估计法的第一步，也是最重要的一步，是构造似然函数。似然函数的基本公式是：
\[L(\theta)=\prod_{i=1}^{n} f\left(x_{i} ; \theta\right)\]
也就是概率密度函数的联乘积。
第二步，取对数，并对参数进行求导，令求导后的式子等于0，解出来的参数就是极大似然估计量。
也就是求解：
\[\frac{\partial}{\partial \theta} \ln L(\theta)=0\]
第三步是验证。把刚才的式子再对参数求一次导数，（即二阶导）若二阶导小于0，就说明求出来的参数的确是最大值，
即验证：
\[\frac{\partial^{2}}{\partial \theta^{2}} \ln L(\theta)<0 \Rightarrow \max\]
下面讲一下对于censored data和truncated data分别如何进行极大似然估计。后两步同上，唯一的区别在于构造似然函数。
Censored data的似然函数
比如说，我们的样本里有n个准确的观测值，和m个删失数据，对这m个删失数据，只知道它们大于y。那么似然函数可以写为：
\[L(\theta)=\left[\prod_{i=1}^{n} f\left(x_{i} ; \theta\right)\right] \times[\mathrm{P}(\mathrm{X}>y)]^{\mathrm{m}}\]
也是密度函数或概率的联乘积。
Truncated data的似然函数
这个我们来看一道例题。
这里我们对于索赔数为0的保单的信息是一无所知的。因此我们在计算各索赔数的保单的概率时，在分母的位置需要扣除索赔数为0的保单。
7
在讲小白鼠的具体例子之前，我们先介绍卡普兰·迈耶生存分析法。
卡普兰·迈耶生存分析法，又称为乘积极限法，是非参数分析法的一种，是一种描述终点事件发生情况的有效手段。Nelson-Aalen模型的原理与之类似，也是生存分析法的一种，用到的是积分风险函数，这里不再介绍，也是CT4第八章的内容。
假定我们观测了n个对象，其中有m个在观测期间发生失效事件（比如我们需要研究的目标事件是作用在小白鼠身上的格列宁直接引起的死亡）。
令t1，t2，…，tk为按顺序观测到的失效事件发生的时刻；
假设在时刻tj发生失效事件的个体数目dj为，则有d1+d2+…+dk=m；
剩余n-m个研究对象处于存活状态的个体是删失的，对这些删失数据，记在时间区间（tj,tj+1）发生失效事件的个数为cj，并设定t0=1、tk+1=无穷，则有c0+c1+…+ck=n-m；
假定失效事件发生的风险在真实发生失效事件时为正，未发生时均为0；
设nj为在时刻tj前存活且能被观测到的数目，定义发生失效事件（这里只研究死亡）的风险为：
\[\hat{\lambda}_{j}=\frac{d_{j}}{n_{j}}(1 \leq j \leq k)\]
离散的死亡风险函数（discrete hazard function of death）为：
\[\lambda_{j}=P\left[T=t_{j} | T \geq t_{j}\right]\]
余命T的分布函数为：
\[F(t)=1-\prod_{t_1\leq t}\left(1-\lambda_{j}\right)\]
余命T的分布函数的估计量为：
\[\hat{F}(t)=1-\prod_{t_{1} \leq t}\left(1-\hat{\lambda}_{j}\right)=1-\prod_{t_{1} \leq t} \frac{n_{j}-d_{j}}{n_{j}}\]
8
回到我们实验室里观测的15只小白鼠的例子。
精算后花园实验室对注射新药格列宁后的15只小白鼠进行了为期30天的死亡调查，希望得到小白鼠注射新药后的生存数据。具体的死亡情况如下图：
我们可以画出时间线：
用上一节的标记方法。本例子中，观测对象总数目n=15。因为注射新药直接死亡的小白鼠数目m=6，分别是Rat 4，Rat 13，Rat 6，Rat 9，Rat 10，Rat 5。
注意，在第17天死亡的小白鼠的死因不是注射新药，而是被其他鼠谋杀，所以这种情况和因为注射新药直接死亡的情况是区别对待的，属于一种随机右删失。
发生失效事件（因为注射新药直接死亡）的时刻分别为t1=3, t2=4, t3=11, t4=21, t5=26。
在对应时刻发生的死亡的小白鼠的数目：d1=1, d2=1, d3=2, d4=1, d5=1。
在对应时刻仍然存活且能被观测到的小白鼠数目：n1=15, n2=14, n3=12, n4=9, n5=6。
删失分别发生在第6天（Rat 7咬坏笼子逃跑了），第17天（Rat 1被其他鼠谋杀），第24天（Rat 8被动物保护主义者突击搜捕后放掉了），第25天（Rat 12被香港记者在采访中放掉了），第30天（剩下的5只老鼠在调查结束后放回了大街上）
所以对应时间区间里删失数据的个数：c0=0, c1=0, c2=1, c3=1, c4=2, c5=5
对应的数据已经找出来，就可以按照上节的公式进行生存分析了。
9
至于那只小白鼠，后来听说它带着血液里的格列宁远渡重洋去了印度，在那里开了一家仿制格列宁的药厂，前一阵子因为一部电影的热映，药价高昂的问题还引起了国内广泛的讨论。
但我最关心的却不是这些。
生老病死的风险，可以通过提前做保险配置，得到一些经济上的补偿，给人更多活下去的勇气以及希望。做一些预备方案，多少是能有点办法的。
但对于这只小白鼠的出逃，我束手无策。
小白鼠离开的第N天。想它。
好久不见，它是不是还喜欢穿粉红色的T恤？它的吃相还是和以前一样喜庆么？
也不知道它长胖了没啊。
10
过去只是我们在观测期间偶然相遇而发生的故事，而这30天里的故事总是绕不开你。
“你我相逢在黑夜的海上
你有你的/我有我的/方向
你记得也好
最好你忘掉
在这交会时互放的光亮”
每一次分别，要是都像处理censored data一样简单，那就好了。
再不回来，我就真把你当censored data了喔。
正在检查 Disqus 能否访问...