随访数据的统计分析方法 完整版

随访数据的统计分析方法

——生存分析与回归模型

第一节、概述：

在医学随访研究中，我们着重收集病人出现某种结果（痊愈、复发、失败或死亡）所经历的时间，以便比较不同干预措施的远期效应的优劣。在这类研究中，有时我们还收集一些有关因素，以分析哪些因素对延长生存时间有利，哪些因素对延长生存时间不利。生存时间的统计分析方法起源于19世纪对寿命表的研究，在第二次世界大战期间，由于对武器的可靠性的要求，使这一分析方法得到了很大的发展，并不断扩展应用的其他研究领域中。近40年来，在医学研究，特别是在临床随访研究中，也引进了生存分析的方法，用来分析病人的随访资料。由于临床研究资料的多样性和复杂性，反过来又进一步推动了生存时间分析技术的发展。到目前为止，生存分析作为统计学的一个分支，已形成了一套完整的体系，包括参数法，非参数法以及回归分析方法等。（一）生存时间（SurvivalTime)：从狭义的角度来说：生存时间是患某病的病人从发病到死亡所经历的时间跨度。开始发病病人死亡t0tjSurvivalTime从广义的角度：从某种起始事件到达某种终点时间所经历的时间跨度。起始事件终点事件t0tjSurvivalTime一、基本概念起始事件终点事件疾病确诊死亡治疗开始痊愈症状缓解疾病复发接触毒物出现毒性反应接触危险因素发病长第1颗牙出现龋齿

生存时间分析时必须要有明确的时间起点和终点以及关于时间的测度单位（如：小时、日、月、年）（二）生存时间数据的类型：1.完全数据某个观察对象具有明确的结局时，该观察对象所提供的关于生存时间的信息是完整的。我们把达到了明确结局的观察对象的生存时间数据称为完全数据（CompleteData）。2.截尾数据由于某种原因未能观察到观察对象的明确的结局，所以不知道该观察对象的确切的生存时间，就象该观察对象的生存时间在未到达规定的终点之前就被截尾了，因此称之为截尾数据（CensoredData）。截尾数据提供了部分关于生存时间的信息，使研究者知道该观察对象至少在已经经历的这个时间长度内没有发生终点事件，其真实的生存时间只能长于我们现在观察到的时间而不会短于这个时间。（符号t+）观察日期t观察起点观察终点确诊死亡截尾病人号123490.1.192.12.3195.12.31t1=2yt2=4+yt3=3yt4=1+y图1完全数据与截尾数据产生截尾现象的原因：1.观察对象的失访，如迁移等原因。2.发生意外，如病人死于他病、动物的意外死亡等。3.研究被终止时观察对象乃未出现终点事件。截尾数据的类型：1.随机截尾(Randomcensoring)可在观察期内的任何时间点上产生截尾。如：3(m)5661014+171920+22+252.右截尾(Rightcensoring)或称定时截尾所有到了某个时点未发生终点的观察对象的生存时间都属于截尾数据。如：4(w)689101011121212+12+12+12+

（三）生存时间资料的特点

（1）效应变量有2个：一是生存时间（天数），二是结局（死亡与否、是否阳性等）。（2）存在截尾数据：如表图中的2号和4号病人未观察到底，不知他们究竟能活多长时间。（3）分布类型复杂：生存时间资料常通过随访获得，因观察时间长且难以控制混杂因素，故其分布常呈偏态，影响因素较多，规律难以估测。

二、统计描述

一般的统计描述属于空间分布范畴，侧重于集中趋势（如平均数、死亡率）和离散趋势（如标准差）。随访资料中时间数据的分布属于时间分布，强调时间过程和截尾数据，故死亡强度必须用概率表示，生存时间不能计算均数，此外还要有将时间与死亡结合在一起的指标。

1、死亡概率（mortalityprobability）：记为q，是指往后一个时段内死亡的可能性大小（概率预测未来的，即将发生的情况）。年死亡概率的计算公式为

q=期内死亡数／期初观察例数若观察期内有截尾，则分母用校正人数，例如，校正人数＝期初人数一（截尾例数／2）死亡率与死亡概率两者的计算和意义都有区别，即①计算：公式中分母不同；②意义：死亡率反映年平均死亡强度，而死亡概率表示往后一年死亡的机会大小。

2、生存概率（survivalprobability）；记为p，是死亡概率的对立面，指往后活满一个时段的可能性大小。年生存概率表示往后再活一年的机会大小，其计算公式为

p＝1-q＝活满t时刻的人数／期初观察例数分子部分是t时刻尚存人数，若年内有截尾，则分母用校正人口数。

t040t130t220t410P0=30/40=0.75p1=20/30=0.667……表1

3、生存率及其标准误（1）生存率（survivalrate）：记为S（t），是指病人经历给定的t个单位时间之后仍存活的概率。生存率是一种习惯叫法，其实质是累积生存概率（cumulativeProbabilityofsurvival〕。若无截尾数据，其公式为

S（t）＝P（T≥t）＝t时刻仍存活的例数／观察总例数其中T为病人的存活时间。上式的定义和意义都很明确，但如果含有截尾数据，则观察总例数将难以确定，故此式一般不能直接应用。

当含有截尾数据时可利用概率乘法原理估计生存率。假定病人在各个时段生存的事件独立，各时段的生存概率为pj，j＝l，2，…，k，则应用概率乘法原理得累积生存概率估计的公式为

S（t）＝P（T≥t）＝p1p2…pk

式中pj可用校正人数估计，可处理截尾数据。上例：3年生存率为

S（3）=10/40=0.250

由式（12-5）求得例12-1的3年生存率为

S（3）＝p1p2p3＝0.750×0.667×0.500＝0.250

（2）生存率的标准误：生存率的标准误有不同的估计方法，其中Greenwood’s法（1926）比较常用，其公式为

（3）生存曲线（survivalcurve）：是指将各个时点的生存率连接在一起的曲线图。）估计的生存率是间断性的，曲线形状分两种类型：①阶梯型：小样本资料用直接法估计的生存曲线（图）；②折线型：大样本资料用频数表法估计的生存曲线（图）。

4、半数生存期及其四分位数间距病人治疗后至死亡所经历的时间跨度即是生存期，但截尾病人的存活时间不是生存期的真值。半数生存期及其四分位数间距的意义同普通的百分位数，但因生存时间资料存在截尾数据，故不能采用普通百分位数的公式计算，一般通过生存率反推时间的方法来估计。

（1）半数生存期（mediansurvivaltime）：又称中数生存期，记为T50，其定义为

T50＝生存率为0.5时所对应的时间它表示有并且只有50％的个体可活这么长时间，它反映生存期的平均水平。如例题12-1，从上列可知，S(2)＝0.500，故其半数生存期是2年，即病人平均可存活2年。当样本生存率恰好有S（tk）＝0.5时，半数生存期直接由定义估计，分两种情况：

若样本生存率中无0.5这个数值，可用插值法估计（不分曲线类型）。例如有两个生存概率分别为：S（T＝90）＝0.4，S（T＝80）＝0.6，则其半数生存期为

（2）四分位数间距：记为Q，表示中间半数病人生存期的分布范围，它反映生存期的离散程度，其定义为

Q＝T25-T75

式中T25和T75，分别是25百分位数和75百分位数，其估计方法同半数生存期，即直接由定义估计或用插值法估计。三、资料要求

生存时间资料协变量的要求同一般统计资料，但对应变量有特殊的要求：

1．死亡例数和比例不能太少这类资料的样本大小主要看死亡例数和比例，而不是总例数，因其信息主要由死亡病例提供，死亡病例少则信息量小，死亡比例小则易出现偏性。

2．截尾原因无偏性例如，老年患者常因不重视随访而失访，由此可能使估计的生存率偏高。为防止截尾偏性，常需对被截尾者的年龄。职业和地区等构成情况进行分析。

3．生存时间尽可能精确因为多数生存分析方法都是在生存时间排序的基础上进行的，即使是小小的合人误差，也可能改变生存时间顺序而影响结果。对于随访资料，生存时间最好精确到天数。

第二节小样本生存率的Kaplan－Meier估计

此法直接用概率乘法原理估计生存率，故又称乘积极限法（Product-limitmethod），它由Kaplan－Meier于1958年提出，故又称Kaplan-Meier法（记为KM法）。这是一种非参数法，主要用于小样本资料的生存率估计。一、生存率估计

例12-2一组病人的存活时间为79，133，185，475，133，238+天，试用Kaplan－Meier法估计生存率。参照表2，计算步骤为：

1.数据列表①将存活天数的数据（1）从小到大排列并列出序号（j），重复数据只列一次（如“133”），第（1）（2）列；②登记各时刻的死亡人数（d）和截尾人数（c）。第（3）、（4）列）

2.求t时刻期初例数n0

本例最后时刻期初人数n5＝1，其它时刻由下往上累计获得，例如，

n03＝n04＋d3＋c3＝2＋1＋0＝3，n02＝n03＋d2＋c2＝3＋2＋0＝5（见第5列）

3.求t时刻死亡概率q＝d／n0。（见第6列）

4.求t时刻生存概率p＝1-q。（见第7列）

5．计算生存率及其标准（见第8，9列）

二、生存曲线

用例12-2资料的生存率所绘制的生存曲线如图所示。此法估计所有死亡时点处的生存率，相邻两个死亡时点之间生存率不变（如本例79～132天的生存率都是0.833），但在右端点处死亡概率即刻改变，生存率即时降低（在133天生存率立即降为0.500），故其曲线有如下特点：

1.曲线阶梯形不能用直线或曲线连接相邻的两个生存率散点。随着死亡时点的增多，曲线的阶梯形就不明显了。

2.曲线左连续每一级台阶的右端为断点（空穴），当前死亡时点处的纵坐标值在下一个台阶。如图中S（133）＝0.5（不是0.833）。

3.曲线尾部不稳定随着时间的增加，观察例数越来越少，误差愈来愈大，曲线尾巴极不稳定。在多组比较时，常发生曲线尾部交叉现象，这很可能是因误差大而出现的一种假象。此时，可适当提前终止日期，使得最后一个死亡时点仍有一定的观察例数（如n＞10）。例如，如果本例将终止时间提前一天，则474天的生存率还是0.333，而475天的生存率就不估计了，新的结果如图，这样曲线尾部的竖线消失，使曲线的规律性更好些。

第三节大样本生存率的寿命表法估计

当样本例数及死亡数足够多时，乘积极限法可按时间分组计算，这就是寿命表法（life-tablemethod）。寿命表法的应用早于乘积极限法，但实质上是乘积极限法的一种近似（频数表法）。寿命表有现时和定群之分，估计生存率用的是定群寿命表法（cohortlifetable）。

一、生存率估计

例12-3某恶性肿瘤随访资料如表3中第1－4列所示，用寿命表法估计生存率。计算方法和步骤为：

1.数据列表①按术后年数从0开始列出时段（t～）及其序号周，②分别列出各时段的死亡和截尾例数。（第1～4列）

2.求年初人数本例存活5年及其以上者24人，各时段年初人数由下往上累计获得，例如，

n05＝n06＋d5＋c5＝24＋9＋0＝33，n04＝n05＋d4＋c4＝33＋18＋1＝52。（第5列）

3.求校正年初人数nc＝n-e/2（第6列）

4.计算死亡概率q＝d／nc（第7列）

5.计算生存概率p＝1-q（第8列）

6.计算生存率及其标准误按（12-5）式和（12-6）式分别计算。（第9，10列）计算生存率及其标准误后，可用正态近似法估计总体生存率的置信区间。如本例手术后2年生存率的95％置信区间为

0.4185±11.96×0.0315，即（0.3567，0.4802）

二、结果分析

1.动态变化：从表3的死亡概率一栏看，前3年各时段死亡危险性逐年增加，而后呈下降趋势；生存概率一栏从反面支持了这一结果。

2.累积情况：从表3的第10列可见，生存率的标准误都很小，说明生存率具有代表性；再看第9列的生存率，半数以上的病人术后活不到2年，结果提示此恶性肿瘤对生命威胁大。

三、生存曲线

用例12-3t年生存率绘制的生存率曲线见图。

Kaplan-Meier法对所有死亡时点估计生存率，其生存率变化是跳跃式的，故曲线呈阶梯形。寿命表法只估计时段右端点的生存率，省略了时段内的生存率估计，但绘图时不能忽略时段内生存率的变化，故曲线不是阶梯形，其特点为：

1.曲线折线形因未知时段内生存率的变化规律，故用直线连接之，组成一条折线。

2.曲线连续可估计任意时点的纵坐标值（生存率）。

3.曲线尾部稳定性较好寿命表法用于大样本，通常最后一个时段仍有一定的观察例数，故曲线尾部的稳定性比Kaplan-Meier法好。

第四节生存曲线比较的假设检验

本节着重讨论两组比较的检验方法，其零假设H0为两总体生存曲线相同，但检验过程一般不估计生存率和半数生存期，而利用死亡数和死亡率函数作统计推断。

一、Logrank此检验的基本思想

Logrank检验的基本思想是实际死亡数与期望死亡数进行比较。例12-4两组淋巴肉瘤患者治疗后复发时间（月数）如表4，对照组为“摘除十放疗”，处理组为“摘除十放疗十化疗”，问在“摘除十放疗”基础上附加“化疗”是否可延长缓解期（相当于寿命资料的生存期）？本例乘积极限法估计的缓解曲线见图，两条曲线分辨度很好。

在H0成立的条件下，此统计量服从自由度为1的χ2分布。对例12-4资料作log-rank检验的计算步骤为：

1.检验假设H0：两总体缓解率曲线相同，H1：两总体缓解率曲线不同；α＝0．05。

2.时间排序将两组非截尾时间混合从小到大排序，整理成多个四格表。约定每个四格表中两组合计的暴露数为Nj，复发数为m1，未复发数为m2，每组的例数分别为n1和n2，计算其中一组的理论数e（本例约定对照组，如果换为处理组也是等价的）。每个四格表的左上格子中的值为aj。

3.求理论数及其方差在H0成立的条件下，每个四格表中的aj均为随机变量（其余3个量本身不是独立变量，只是随aj变化而已），其理论数与方差为:

例：t=2(月):结局对照处理合计复发2(a)13(m1)缓解121628(m2)合计14(n1)17（n2)31(N)Na=31*2=62Ne=n*D=14*342N2v=312*0.693=665.973

4.求统计量当这些四格表足够多时，将所有aj合并得实际数合计∑aj，理论数ej合计∑ej，方差vj合计∑vj，可用正态近似法作χ2检验。Log-rank检验统计量为

得∑aj＝7，∑ej＝3.820，∑vj＝2.229，代入（12-11）式得

5.确定P值并作出结论据自由度为1的χ2分布求得P＝0.0332，按α＝0.05水准拒绝H0，接受H1，故可认为附加放疗可使患者尽早缓解。

二、Breslow检验

如果随机变量aj以四格表的例数Nj为权重，便可获得Breslow统计量为

在H0成立的条件下，此统计量亦服从自由度为1的χ2分布。对例12-4资料作Breslow检验的计算步骤为：

1~3步与上面介绍的log-rank检验相同，对每一时刻用Nj加权求的Na、Ne及N2v。

4.求统计量得：∑Njaj＝176，∑Njej＝98，∑Nj2vj＝1473.178，代入（12-12）式得

5.确定P值并作出结论据自由度为1的χ2分布求得P＝0.0421，按α＝0.05水准拒绝H0，接受H1，故可认为附加放疗可使患者尽早缓解。

三、组间差别大小的度量

上述Log-Rank与Breslow检验不计算描述性统计指标，当检验结果有统计学意义时，可通过以下几个途径来考察效果的好坏和差别的大小。

1．生存曲线图（或表）目测判断如从图可知，处理组缓解率高于对照组。

2．半数生存期比较本例Md对照＝33月，Md处理＞69月（因S（69）＝0.816，往后就截尾了），说明处理组的缓解期比对照组长。前两类指标属于生存情况的描述指标，在本章第一节已作介绍。

3．相对危险度RR（RelativeRisk）比较相对危险度是因素效应指标，反映因素作用的相对大小。在作多组比较时，第i组与第j组的相对危险度为

式中Ai、Ei及Aj、Ej分别为第i组与第j组的全部实际死亡数与期望死亡数。例中已获得两组淋巴肉瘤治疗后复发率差别有统计学意义的结论，求处理组（1组）相对于对照组（0组）的相对危险度。

解已知A0＝7，E0＝3.820，A1＝3，E1＝（A0＋A1）－E0＝（7+3）-3.820＝6.180，代入得：

说明处理组危险度是对照组的26.5％，“处理”是保护因子，即在“摘除十放疗”基础上附加“化疗”可望提高缓解率（1／0.265＝3.8倍）。

在生存曲线比较分析中要注意以下几点：（1）方法选择本例上述两种检验的结论一致，相差甚微，但在应用上是有差别的。Log-rank法将所有四格表一视同仁，因四格表例数是逐渐减少的，所以它实际上相对重视了远期效应；而Breslow法以四格表的例数为权重，则相对重视了近期效应。要根据对近期效应和远期效应的重视程度来选择方法。实际应用中常同时计算这两种统计量，当两者结论一致时，可认为近期与远期疗效都有差别（或都无差别）；若Breslow法拒绝H0而log-rank法不拒绝H0，则说明近期疗效有差别而远期疗效无差别，反之亦然。

（2）应用条件除了生存资料的基本要求之外，还要求各组生存曲线不能交叉。若出现这种交叉，则提示可能存在混杂因素，应采用多因素方法来校正混杂作用或分段作统计分析。tS(t)第五节、Cox比例风险回归简介英国统计学家Cox于1972年提出这一模型，用来分析常有伴随变量的生存时间资料。由于这一模型是以生存时间的长短顺序作为分析基础，对生存时间的分布形式没有事先的假定，因此适应性较强，在生存时间资料的分析中的应用比较广。th(t)h0h1h0h1th(t)非比例风险比例风险（1）模型结构；设有n名病人，第i名病人的生存时间为ti，同时该病人具有一组伴随变量xi1，xi2…xip该病人生存到时间ti的风险函数hi(t）是其基础风险函数ho(t)与相应伴随变量的函数的乘积，其数学表达式为：hi（t）=ho（t）f（β1xi1+….+βpxip）Cox建议伴随变量的函数为指数形式，故Cox比例风险回归可写为：hi（t）=ho（t）exp（β1xi1+….+βpxip）

式中的ho（t）是当所有伴随变量xi1，xi2，…xip都处于0或标准状态下的风险函数，是一个不确定的值，β1，β2…，βp为回归系数，须用实际资料来估计。将ho（t）移至等式左边并去自然对数得：等式左边的部分为相对风险度的自然对数值，等式右边部分为伴随变量与相应回归系术的线性组合。βj（j＝1，2，…，p）的实际意义是：在p-1个伴随变量为一定值时，当伴随变量xj每改变一个测定单位时所引起的相对风险度自然对数值的改变量。各自变量xj的回归系数βj与危险度间呈指数函数关系，这是Cox模型的一个假定条件。当βj＝0时，说明xj对危险度不起作用；βj为正值时xj为危险因子，增大了危险度；βj为负值时xj＝1与xj＝0的两个危险度相比，则得到一个与h0(t)无关的比值，称为相对危险度。（2）实例：表是26例Ⅲ期浆液性卵巢上皮癌病人经手术后的生存时间（月）资料，同时记录了这些病人的年龄（x1＝年龄/10岁）、细胞分化度（x2）、淋巴细胞浸润度（x3）及手术残留灶（x4，单位：cm）。不变量为1时表示死亡，为0时表示生存。按α=0.05水平，用Cox比例风险回归拟合得到对生存时间有影响的因素为：①手术残留灶x4，其回归系数估计值为表示手术残留灶每增加1cm，则死亡的相对风险度上升exp（0.9446）＝2.57（倍）；②细胞分化度x2，其回归系数估计值为表示细胞分化度每增加一个等级时，相对风险度上升exp（0.8205）＝2.27（倍）。而年龄x1与淋巴细胞浸润x3未引入方程。本例Cox模型的参数估计值及相应的方差见表7.10（3）COX模型的参数估计；由于该模型中h0（t）项属于半参数模型，故在构成偏似然函数后，应用最大似然法来估计βj。对预后因素的选择常用后退法，即先把全部预后因素引入模型，然后再通过假设检验将未达α水平的预后因素剔除。常用的假设检验方法有梯度检验、Wald检验（即u检验）与似然比检验，在大样本检验中三者的结论是一致的。

（4）Cox模型在临床试验中的应用：Cox模型适用处理单因素或多因素影响下的时间－效应数据，这种数据除包括每个研究对象的有关因素外，还有生存时间，或结局未发生的截尾数据。在此有必要说明：仅有截尾数据是不能进行时间－效应分析的，截尾数据必须结合有结局发生的完全数据。截尾数据所占比重太多，建立的模型也不稳定。在临床试验的定群研究和随访研究中，常用Cox模型作为分析的重要手段，尤其是在研究干预措施的远期效果，在设计方面比其他模型更显得灵活与适用。如对随访迟早不一、随访时间长短不一、夫访的部分信息等方面都容易实施。然而此法的计算和繁复的公式对初学者有一定难度，但现在随着电脑知识的普及和统计软件的开发，这些难点都将迎刃而解。表140例癌症病人生存分析基本指标计算表（无截尾）术后年初年内死年内截年平均年死亡数年死亡概率年生存例数t+1年生存率年数例数亡例数尾例数例数t～n0dcn=n0-d/2m=d/nq=d/n0p=(n0-d)/n0S(t+1)=(n0-d)/N(1