第十一章 生存分析在医学研究中有时候需要对病人的生存情况加以评价，例如：肺癌病人术后生存情况或者是两种不同方案对肺癌病人治疗后的比较。从理论上说，在完全收集到所有病人因肺癌死亡的确切时间后，可以用术后生存时间这个指标来描述病人的生存状况，并对不同的组别间加以比较。但在实际随访过程中，由于失访、意外死亡等原因，部分研究对象不能随访到确切的存活时间。所以，无法以生存时间直接进行统计分析，而需要采取一些特殊的统计方法，接下来将对此进行详细介绍。
基本概念
1．1 生存时间、截尾例11.1：某医师从2002年1月1日起对６名肺癌患者进行跟踪观察，记录其结局，随访记录见表11.1：
表11.1 6例肺癌患者随访记录研究号
(1)
姓名
(2)
术后开始随访时间
(3)
终止随访时间
(4)
结局
(5)
生存天数
(6)
1
赵××
00-01-16
00-09-01
失去联系
229+
2
钱××
00-01-26
01-03-18
局部复发死亡
417
3
孙××
00-02-05
01-12-25
车祸死亡
689+
4
李××
00-02-19
02-10-01
研究终止
955+
5
王××
00-03-07
00-11-11
局部复发死亡
249
6
刘××
00-03-20
02-01-24
远处转移死亡
675
生存时间（survival time）是指从某起点事件(start point event)开始到被观测对象出现终点事件（endpoint event）所经历的时间。如从疾病“确诊”到“死亡”；从“治疗结束”到“复发”；某电子设备从“开始使用”到“出现故障”，等。由此可见，此处的“生存”是一个广义的概念。根据研究目的不同，可以有不同的“起点时间”和“终点时间”。生存时间的单位可以用年，月，周，日，甚至时，分，秒等表示。上表中的(6)即为生存时间，有2种类型：
（1）．完全数据（complete data）：是指被观测对象从观察起点至出现终点事件所经历的时间。在例11.1中，如果终点事件被设为死于肺癌，那么417、249、675为完全数据。
（2）．截尾数据（censored data）：是指在出现终点事件前，被观测对象的观测过程终止了。由于被观测对象所提供的信息是“不完全的”，只知道他们的生存时间超过了截尾时间。例11.1中，229+、689+、955+为截尾数据。
起点事件（start point event）：指研究者根据研究目的开始关心某一事件的起点，如上面所说的“确诊”、“治疗结束”、“开始使用”等等。
终点事件（endpoint event）：指研究者根据研究目的所关心的某一事件，如上面所说的“死亡”、“复发”、“出现故障”等等。需要注意一点，不同的研究目的有不同的终点事件，如果研究的是肿瘤的局部复发情况，那么死于肿瘤远处转移只能算做截尾，而不是终点事件。
截尾(censor)，又称终检、删失，主要有3种情况：①失访：指失去联系，如信访未回信，上门不见人，电话采访不答理或搬迁未留新地址等原因；②退出：指退出研究，如因其他非此次研究疾病死亡，临时改变治疗方案而中途退出研究；③终止：指研究时限已到而终止观察。
1．2 死亡概率、生存概率
死亡概率（mortality probability）：记为q，是指在某段时间开始时存活的个体在该段时间内死亡的可能性大小。若无截尾数据，死亡概率的估计公式为：
(11.1)
生存概率（survival probability）：记为p，是指在某段时间开始时存活的个体至该时间结束时仍存活的可能性大小。生存概率的公式为：
(11.2)
由于生存分析中常存在截尾值，假定失访等截尾事件在观察时段的各个时间点等机会发生的，故分母改用校正观察例数
1．3 生存率及标准误、生存曲线生存率（survival rate）：用S(tk)表示，是指病人经历tk个单位时间后仍存活的概率。若无截尾数据，则
(11.3)
其中t为病人的存活时间，但如果资料中含有截尾数据，分母必须按时段进行校正，此公式则不再适用，此时生存率的计算公式应为
(11.4)
其中、、…、表示不同时间段的生存概率，可以看出，生存率是多个时段生存概率的累积，故生存率又称为累积生存概率（cumulative probability of survival）。
生存率标准误的计算公式为
(11.5)
例11.2：某医院对100名恶性肿瘤术后病人进行随访所获得的资料如表11.2：
表11.2 某恶性肿瘤100例术后随访资料术后年数
t～
(1)
期初观察例数
n0
(2)
期内死亡例数
d
(3)
期内截尾例数
c
(4)
校正期初例数
nc=n0-c/2
(5)
死亡概率
q=d/n
(6)
生存概率
p=1-q
(7)
t+1年生存率
S(t+1)
(8)
生存率标准误
SE(S(t+1))
(9)
0～
100
10
0
100
0.1000
0.9000
0.9000
0.0300
1～
90
15
0
90
0.1667
0.8333
0.7500
0.0433
2～
75
16
8
71
0.2254
0.7746
0.5810
0.0501
3～
51
12
4
49
0.2449
0.7551
0.4387
0.0520
4～
35
14
7
31.5
0.4444
0.5556
0.2437
0.0484
5～
14
7
7
10.5
0.6667
0.3333
0.0812
0.0390
（1）为组段，（2）为各组段期初观察例数，（3）、（4）分别为各组段内的死于该恶性肿瘤例数和截尾例数，（5）为校正期初例数，（6）、（7）分别为各组段的死亡概率和生存概率，计算公式分别为(6)=(3)/(5)、(7)=1－(6)，（8）为t+1年生存率，（9）为t+1年生存率的标准误。
在本例中，q1=10/100，p1=1-q1=90/100，q2=15/90，p2=1-q2=75/90。
由式（11.4）求得2年生存率为S(2)= p1p2=(90/100)×(75/90) =75/100=0.75
由式（11.3）求得2年生存率为S(2)=(100－10－15)／100＝75/100=0.75
两式求得结果相同，但在继续计算S(3)、S(4)…时，由于式（11.3）不能对截尾数据进行处理，故不再适用，只能以（11.4）进行计算。
生存曲线（survival curve）：以时间为横轴，生存率为纵轴，将各个时点的生存率连接在一起的曲线图。
1．4 中位生存期中位生存期(Median Survival Time)：又称为半数生存期，即当累积生存率为0.5时所对应的生存时间，表示有且只有50%的个体可以活过这个时间。要注意中位生存期通常不等于生存时间的中位数（除非在这个时间点之前没有删失值存在）。
第二节 生存曲线估计
对于生存曲线的估计，我们介绍常用的2种方法：
1．Kaplan-meier法，又称乘积极限法（Product-Limit Method，简称PL法），由Kaplan-Meier在1958年提出，适用于样本量较小，难以将生存时间按组段划分，此时是利用tk时刻之前各时点上生存概率的乘积来估计在时刻tk的生存率，不需要对被估计的资料分布作任何假设。
例11.3：一组病人的生存时间（日）如下，用Kaplan-Meier法估计其生存曲线（+代表截尾）。
3,5,5,6+,8,16+,22,30,47+,71
为便于说明，将计算过程和结果列在表11.3，其步骤为：
（1）．将所有生存时间按从小到大排列（t），包括完全和截尾生存时间。
（2）．列出各期初暴露病例数（n），它是指在t时刻前仍存活的病例数。
（3）．将各期内死亡例数（d）和截尾例数（c）分别写在第（3）、（4）列。
（4）．计算各期的死亡概率q,q=d/n，截尾数不参加计算，结果见第（5）列。例如生存时间为5天时，q=2/9=0.2222。
（5）．计算各期的生存概率p,p=1-q，例如生存时间为5天时，p=1-0.2222=0.7778，计算结果见第（6）列。
（6）．计算各t时刻的生存率S(tk)。计算tk时刻生存率时可以用小于和等于tk时刻的各时点生存概率的乘积得到，计算结果见第（7）列。
表11.3 Kaplan-Meier法估计生存率计算表存活时间（天）
期初暴露病例数
期内死亡例数
期内截尾例数
死亡概率
生存概率
生存率
生存率标准误
t
n
d
c
q=d/n
p=1-q
S(tk)
SE(S(tk))
(1)
(2)
(3)
(4)
(5)
(6)
(7)
(8)
3
10
1
0
0.1000
0.9000
0.9000
0.0949
5
9
2
1
0.2222
0.7778
0.7000
0.1449
6
7
0
1
0.0000
1.0000
0.7000
0.1449
8
6
1
1
0.1667
0.8333
0.5833
0.1610
16
5
0
1
0.0000
1.0000
0.5833
0.1610
22
4
1
0
0.2500
0.7500
0.4375
0.1747
30
3
1
1
0.3333
0.6667
0.2917
0.1666
47
2
0
1
0.0000
1.0000
0.2917
0.1666
71
1
1
0
1.0000
0.0000
0.0000
.
用Stata软件计算时，数据库结构如下：
time
dead
1
3
1
2
5
1
3
5
1
4
6
0
5
8
1
6
16
0
7
22
1
8
30
1
9
47
0
10
71
1
Stata命令如下：
stset time,failure(dead=1)
将该数据库定义为生存资料数据库
sts list
输出各时点生存率
输出结果：
Beg,Net Survivor Std.
Time Total Fail Lost Function Error [95% Conf,Int.]
-------------------------------------------------------------------------------
3 10 1 0 0.9000 0.0949 0.4730 0.9853
5 9 2 0 0.7000 0.1449 0.3287 0.8919
6 7 0 1 0.7000 0.1449 0.3287 0.8919
8 6 1 0 0.5833 0.1610 0.2298 0.8207
16 5 0 1 0.5833 0.1610 0.2298 0.8207
22 4 1 0 0.4375 0.1747 0.1187 0.7256
30 3 1 0 0.2917 0.1666 0.0471 0.6085
47 2 0 1 0.2917 0.1666 0.0471 0.6085
71 1 1 0 0.0000,,,
-------------------------------------------------------------------------------
输出结果依次为时点、各时点的期初例数、死亡例数、截尾例数、生存函数（生存率）及其标准误、生存率的95%可信区间。如果需要输出生存函数曲线图，接着输入命令：
sts graph,title（Kaplan-Meier法对10名病人的生存率估计）
输出生存曲线
输出结果：
图11.1 Kaplan-Meier法对10名病人的生存率估计其中title()选项要求输出图像的标注。
2．寿命表法（Life-Table Method，简称LT法）是通过计数落入时间区间［tk-1,tk）内的失效和截尾的观察例数来估计该区间上的死亡概率，然后用该区间及其之前各区间上的生存概率之积来估计S(tk)。当样本量较大或者无法准确得知研究结果出现的时间时，可以将各研究对象的生存时间按某个时间段（年、月等）进行分组计算其生存率。
对例11.2用寿命表法计算生存率，具体计算步骤前面已有叙述，这里主要介绍stata的用法。
数据库结构如下：
time
num
dead
1
0
10
1
2
1
15
1
┆
┆
┆
┆
6
5
7
1
7
0
0
0
┆
┆
┆
┆
11
4
7
0
12
5
7
0
此时的数据库中有3个变量，time为生存时间（年），time=n表示该组患者在第n年初随访到了但是在n+1年初没有随访到。dead是观察结局，dead=1表示死亡，dead=0表示截尾。num为各组频数。
Stata程序如下：
use d:\data\ltable1
打开ltable1数据文件
ltable time dead [weight=num]
制作寿命表
输出结果：
Beg,Std.
Interval Total Deaths Lost Survival Error [95% Conf,Int.]
-------------------------------------------------------------------------------
0 1 100 10 0 0.9000 0.0300 0.8221 0.9449
1 2 90 15 0 0.7500 0.0433 0.6529 0.8236
2 3 75 16 8 0.5810 0.0501 0.4765 0.6718
3 4 51 12 4 0.4387 0.0520 0.3354 0.5371
4 5 35 14 7 0.2437 0.0484 0.1557 0.3425
5 6 14 7 7 0.0812 0.0390 0.0260 0.1779
-------------------------------------------------------------------------------
输出结果依次为各段生存时间起点与终点、期初人数、期内死亡人数、截尾值人数、生存率、标准误及其相应的95%可信区间。要输出生存函数曲线图，只须在ltable命令后加入graph选项：
ltable time dead [weight=num],graph notab noconf title (100例恶性肿瘤患者术后生存率情况)
图11.2 100例恶性肿瘤患者术后生存曲线本命令中，graph要求输出生存函数曲线；notab要求不输出寿命表；noconf要求在生存函数曲线中不显示95%可信区间，title（）则是给图形加上标注。
如果获得资料中每一例患者都有自己的生存时间和随访结局，并且样本量较大时，也可以用寿命表法进行计算，此时的stata数据库结构为：
time
dead
┆
┆
┆
stata命令为：
ltable time dead,interval(365)
以365天（1年）为组段制作寿命表
得到结果是一致的。
3．两种方法的比较：
①．寿命表法适用于大样本或无法准确得知研究结果出现时间的资料，Kaplan-Meier法主要用于小样本，也可以用于大样本。
②．寿命表法是按照指定的时段来分段，估计的是时间区间右端点上的生存率；Kaplan-Meier法是根据死亡时点分段，逐个估计死亡时点的生存率。
③．寿命表法不能确切得知死亡时间，假定每个时间段中的“死亡”是呈均匀分布，生存率为线性变化，故简单化以直线相连接；Kaplan-Meier法其生存曲线是左连续的阶梯型曲线，间断点的纵坐标在下一阶处，当样本量较大及死亡时点较多时，阶梯形就不明显了。
4．中位生存期的估计
寿命表法由于默认组段内生存率的变化是均匀的，因此可以直接在生存曲线上进行内插（图11.3a），如例11.2，3年生存率为0.5810，4年生存率为0.4387，当生存率为0.5000时，中位生存时间（年）为=3.57。
而在Kaplan-Meier法中，由于估计的是时点生存率，生存曲线是阶梯形的，对中位生存期的估计有两种观点：一种观点认为中位生存期为生存率降到0.5或以下的首个生存时间，另一观点认为需要先将生存率为0.5两侧左端点连线再进行内插。例11.3在两种观点下中位生存期分别为22天和16天（图11.3b）。stata软件默认的中位生存期估计为前一种观点，可通过stsum命令得到实现。
图11.3a 寿命表法中位生存期估计
图11.3b Kaplan-Meier法中位生存期估计
第三节 生存曲线的比较这一节将介绍两组生存资料比较最常用的时序检验（log-rank test），无效假设H0为两条总体生存曲线相同。如果H0成立，两组生存资料来自同一总体，用两组合并的资料估计一条生存函数曲线，该生存函数在各时段中所计算的理论死亡人数与实际死亡人数相差不会太大，否则拒绝无效假设，接受备择假设，认为各组总体生存曲线不同或不全相同。
例11.4：某医生收集23例晚期肺癌患者在接受化疗后的生存时间t（月），按接受治疗方案的不同划分为2组（1为常规方案，2为新方案），问不同的治疗方案对其生存时间长短的影响有无显著性差异。
常规方案组：
3
3+
4
6
6+
7+
8
9+
14
19
新方案组：
5
9
9+
10+
11
11
12
13+
17
18+
19+
20
32+
计算步骤：
（1）建立假设
H0：两总体生存函数曲线相同
H1：两总体生存函数曲线不同
（2）按时间排序将两组未截尾的完全生存时间从小到大混合排序，见表11.4第（2）列表11.4 两组肺癌患者资料分析表序号
生存时间
死亡数
期初观察数
理论死亡数
实际数－死亡数
j
tj
d1j
d2j
n1j
n2j
e1j
e2j
d1j-e1j
d2j-e2j
(1)
(2)
(3)
(4)
(5)
(6)
(7)
(8)
(9)
(10)
1
3
1
0
10
13
0.435
0.565
0.565
-0.565
2
4
1
0
8
13
0.381
0.619
0.619
-0.619
3
5
0
1
7
13
0.350
0.650
-0.350
0.350
4
6
1
0
7
12
0.368
0.632
0.632
-0.632
5
8
1
0
4
12
0.250
0.750
0.750
-0.750
6
9
0
1
3
12
0.200
0.800
-0.200
0.200
7
11
0
2
2
9
0.364
1.636
-0.364
0.364
8
12
0
1
2
7
0.222
0.778
-0.222
0.222
9
14
1
0
2
5
0.286
0.714
0.714
-0.714
10
17
0
1
1
5
0.167
0.833
-0.167
0.167
11
19
1
0
1
3
0.250
0.750
0.750
-0.750
12
20
0
1
0
2
0.000
1.000
0.000
0.000
合计
6
7
3.272
9.728
2.728
-2.728
（3）将不同生存时间的死亡数按组别归入（3）、（4），期初观察数按组别归入（5）、（6），计算时，先将两组每一生存时间tj的资料列成一个2×2表的形式，见表(11.5)，
表11.5 2×2表形式
死亡
生存
合 计
常规方案组
d1j
n1j-d1j
n1j
新方案组
d2j
n2j-d2j
n2j
合 计
Dj
Sj
Nj
然后按照四格表χ2 检验的方法计算理论死亡数。即
eij = ,i=1,2 (11.6)
将不同生存时间的理论死亡数计算结果列在（7）（8）。常规方案组的理论死亡数合计是3.272，新方案组的理论死亡数合计是9.728。全部实际死亡数与理论死亡数之差为：
O1 - E1 = 2.728
O2 - E2 =-2.728
由结果可看出，两组的实际数与理论数之差值是一样的，只是符号相反。因此计算时，选择任一组计算的结果就可以。
对于两组生存率进行Log Rank 检验，其统计量服从χ2 分布，计算公式为：
,＝1 (11.7)
其中
= ,i =1,2 (11.8)
本例中，=2.260，故
χ2 === 3.29
式（11.7）的近似公式为
（11.9）
本例中，3.04，结果略有不同。
下面用Stata软件进行计算：
数据库结构：
time
treat
dead
1
3
常规方案组
1
2
3
常规方案组
0
┆
┆
┆
┆
9
9
常规方案组
0
10
9
新方案组
1
┆
┆
┆
┆
22
20
新方案组
1
23
32
新方案组
0
Stata命令：
use d:\data\km2
打开km2数据库
stset time,failure(dead)
定义为生存数据库
sts list,by(treat)
按组别输出不同时点生存率
sts test treat,logrank
对两总体生存率曲线进行logrank检验
sts graph,by(treat)
按组别输出生存曲线
输出结果：
failure _d,dead
analysis time _t,time
Beg,Net Survivor Std.
Time Total Fail Lost Function Error [95% Conf,Int.]
-----------------------------------------------------------------------------
常规方案组
3 10 1 1 0.9000 0.0949 0.4730 0.9853
4 8 1 0 0.7875 0.1340 0.3809 0.9426
6 7 1 1 0.6750 0.1551 0.2906 0.8825
7 5 0 1 0.6750 0.1551 0.2906 0.8825
8 4 1 0 0.5063 0.1868 0.1396 0.7903
9 3 0 1 0.5063 0.1868 0.1396 0.7903
14 2 1 0 0.2531 0.2019 0.0138 0.6438
19 1 1 0 0.0000,,,
新方案组
5 13 1 0 0.9231 0.0739 0.5664 0.9888
9 12 1 1 0.8462 0.1001 0.5122 0.9591
10 10 0 1 0.8462 0.1001 0.5122 0.9591
11 9 2 0 0.6581 0.1407 0.3200 0.8576
12 7 1 0 0.5641 0.1488 0.2436 0.7928
13 6 0 1 0.5641 0.1488 0.2436 0.7928
17 5 1 0 0.4513 0.1560 0.1549 0.7121
18 4 0 1 0.4513 0.1560 0.1549 0.7121
19 3 0 1 0.4513 0.1560 0.1549 0.7121
20 2 1 0 0.2256 0.1776 0.0151 0.5896
32 1 0 1 0.2256 0.1776 0.0151 0.5896
-----------------------------------------------------------------------------
failure _d,dead
analysis time _t,time
Log-rank test for equality of survivor functions
| Events Events
treat | observed expected
------------+-------------------------
常规方案组 | 6 3.27
新方案组 | 7 9.73
------------+-------------------------
Total | 13 13.00
chi2(1) = 3.29
Pr>chi2 = 0.0696
图11.4 两组肺癌患者的生存曲线比较
从输出结果可以看出，logrank检验结果为=3.29，=1，P=0.0696，按水平，没有理由拒绝H0假设，认为两组的总体生存曲线相同，在本例中，即在晚期肺癌中，没有理由认为新化疗方案和常规化疗方案的疗效不同。
需要注意的是，在做Log-Rank检验的时候，除了生存资料的基本要求之外，还要求各组生存曲线不能交叉。因为这种交叉提示可能有混杂因素，应该采用分段分析或者采用多因素方法来分析。
第四节 指数模型和Weibull模型设生存时间T，如果生存率S(t)满足下列模型：
其中回归系数，协变量，
因此
称生存时间T服从指数分布模型，显然越大，也越大，因此生存率就越小。故称为指数回归模型的风险函数。称为指数回归模型的累积风险函数。
如果生存率S(t)满足下列模型：
则称生存时间T服从Weibull回归模型，显然当(＝1时，Weibull回归模型就是指数回归模型，因此指数回归模型是Weibull回归模型的特例。Weibull回归模型的累积风险函数为，相应的风险函数为。
以上可以归纳为累积风险函数为－ln(S(t))，并且Weibull模型的累积风险函数为－ln(S(t))＝表示累积风险以t的幂函数增长。特别指数分布模型的累积风险以t直线增长。
以下举一个例子说明这类模型的应用。
例11.5：某医生在研究2种治疗方案对肝癌术后患者的疗效，收集了60名患者的生存资料，随机分为两组，分别采用2种不同的治疗方案，同时收集了患者的年龄情况，生存时间的单位为年，死亡的结局定义为1，生存的结局定义为0，资料如下：
方案A(drug=0)
方案B(drug=1)
年龄
结局
生存时间
年龄
结局
生存时间
71
1
0.5
54
0
26.6
70
1
0.1
69
1
0.1
56
1
2.1
66
1
0.6
52
1
0.2
68
1
1.6
67
1
0.4
70
1
0.7
60
1
1.1
53
1
3.1
74
1
0.1
68
1
3.1
53
1
1.3
51
1
2.4
67
1
0.1
66
1
1
53
1
1.8
76
1
1.7
79
1
0.1
57
1
1.1
56
1
0.6
51
1
0.8
57
1
3.6
62
1
7.2
54
1
0.3
76
1
2.3
72
1
0.3
65
1
4.4
61
0
3.5
75
1
2.9
61
1
0.8
56
1
2.2
79
1
0.4
66
1
1.1
67
1
0.2
57
1
1.2
73
1
0.7
78
1
0.1
52
1
0.3
58
1
3.4
74
1
1.2
53
0
2.9
57
1
0.8
62
1
2.4
63
1
0.3
71
1
4
53
1
0.6
76
1
0.8
56
1
2.3
63
1
0.2
51
1
2
62
1
1.5
55
1
2
76
1
0.3
58
1
0.4
51
1
1.9
68
1
0.1
70
1
0.6
Stata数据格式如下
drug
age
d
t
0
71
1
0.5
0
70
1
0.1
0
56
1
2.1
0
52
1
0.2
0
67
1
0.4
0
60
1
1.1
0
74
1
0.1
0
53
1
1.3
0
67
1
0.1
0
53
1
1.8
0
79
1
0.1
0
56
1
0.6
0
57
1
3.6
0
54
1
0.3
0
72
1
0.3
0
61
0
3.5
0
61
1
0.8
0
79
1
0.4
0
67
1
0.2
0
73
1
0.7
0
52
1
0.3
0
74
1
1.2
0
57
1
0.8
0
63
1
0.3
0
53
1
0.6
0
56
1
2.3
0
51
1
2
0
55
1
2
0
58
1
0.4
0
68
1
0.1
1
54
0
26.6
1
69
1
0.1
1
66
1
0.6
1
68
1
1.6
1
70
1
0.7
1
53
1
3.1
1
68
1
3.1
1
51
1
2.4
1
66
1
1
1
76
1
1.7
1
57
1
1.1
1
51
1
0.8
1
62
1
7.2
1
76
1
2.3
1
65
1
4.4
1
75
1
2.9
1
56
1
2.2
1
66
1
1.1
1
57
1
1.2
1
78
1
0.1
1
58
1
3.4
1
53
0
2.9
1
62
1
2.4
1
71
1
4
1
76
1
0.8
1
63
1
0.2
1
62
1
1.5
1
76
1
0.3
1
51
1
1.9
1
70
1
0.6
操作如下：
定义生存分析数据集
stset t,f(d)
指数回归分析命令
streg drug age,d(e) nohr
No,of subjects = 60 Number of obs = 60
No,of failures = 57
Time at risk = 110.4000007
LR chi2(2) = 33.64
Log likelihood = -91.361189 Prob > chi2 = 0.0000
------------------------------------------------------------------------------
_t | Coef,Std,Err,z P>|z| [95% Conf,Interval]
-------------+----------------------------------------------------------------
drug | -1.230369,2670223 -4.61 0.000 -1.753723 -.7070145
age |,0672908,0157282 4.28 0.000,0364642,0981175
_cons | -4.043394 1.001847 -4.04 0.000 -6.006979 -2.07981
------------------------------------------------------------------------------
A方案drug=1与B方案drug=0的风险函数比(hazard ratio)
，
同理年龄增加一岁的风险函数比
也可以利用Stata直接计算得到：streg drug age,d(e)
No,of subjects = 60 Number of obs = 60
No,of failures = 57
Time at risk = 110.4000007
LR chi2(2) = 33.64
Log likelihood = -91.361189 Prob > chi2 = 0.0000
------------------------------------------------------------------------------
_t | Haz,Ratio Std,Err,z P>|z| [95% Conf,Interval]
-------------+----------------------------------------------------------------
drug |,2921849,0780199 -4.61 0.000,1731282,4931142
age | 1.069606,0168229 4.28 0.000 1.037137 1.103092
------------------------------------------------------------------------------
由于指数回归模型是Weibull回归模型的特例，所以通常用Weibull回归模型要优于指数回归模型。Stata命令如下
streg drug age,d(w) nohr
相应的结果如下
Weibull regression -- log relative-hazard form
No,of subjects = 60 Number of obs = 60
No,of failures = 57
Time at risk = 110.4000007
LR chi2(2) = 26.04
Log likelihood = -91.355743 Prob > chi2 = 0.0000
------------------------------------------------------------------------------
_t | Coef,Std,Err,z P>|z| [95% Conf,Interval]
-------------+----------------------------------------------------------------
drug | -1.24424,2984532 -4.17 0.000 -1.829198 -.6592825
age |,0680308,0172641 3.94 0.000,0341938,1018678
_cons | -4.091576 1.104093 -3.71 0.000 -6.255558 -1.927593
-------------+----------------------------------------------------------------
/ln_p |,0105058,1003806 0.10 0.917 -.1862365,2072482
-------------+----------------------------------------------------------------
p | 1.010561,1014407,8300773 1.230288
1/p |,9895491,0993315,8128179 1.204707
------------------------------------------------------------------------------
本例资料用Weibull回归模型进行参数估计与指数模型的参数估计非常接近，说明本例资料服从指数分布的回归模型。
指数分布的回归模型与Weibull分布的回归模型均属于生存分析中的比例风险参数模型，即：协变量的参数表达式与生存时间呈比例，更一般的比例风险模型为
其中回归系数，协变量，
因此
当X的各个分量全为0时，，故称h0(t)为基准风险函数(Baseline Hazard function)。
生存分析的模型可以分为参数模型和半参数模型，指数分布的回归模型和Weibull分布的回归模型都属于参数模型，这些模型可以直接估计生存率，但是要求资料服从对应的分布。这对实际研究的要求很高，往往难以确认。对于研究生存问题的危险因素而言，可以考虑下列的半参数模型：Cox模型。这种模型对资料的要求大大降低了。
第五节 Cox回归
由于log-rank仅能分析一个因素，因此，在同时分析2个或2个以上因素对生存时间影响的时候就显得无能为力了，这时我们就需要通过Cox比例风险模型来解决这些问题。以例11.5加以说明：
例11.5：某医生在研究2种治疗方案对晚期肺癌的影响，收集了34名患者的生存资料，同时收集了患者的年龄因素，资料如下：
方案1治疗组(N=20)
方案2治疗组(N=14)
生存月数
年龄
生存月数
年龄
1
66
6
72
1
70
6+
70
2
64
7
63
3
57
9+
61
4
61
10+
54
4
72
11+
66
5
68
13
67
5
63
15+
55
8
61
16
72
8+
63
19+
55
8
57
20+
60
8
54
22
63
11
55
23
52
11
60
32+
57
12
54
12
67
15
56
17
54
22
62
23
57
希望回答以下问题：
①在同样的年龄情况下，不同的治疗方案对生存时间有无影响？
②在同一治疗方案的情况下，年龄对生存时间有无影响？
③治疗方案和年龄有无交互作用？即在不同的年龄段中，不同的治疗方案对生存期的影响是否不同？
1.模型的结构
在Cox比例风险模型中，假定有一个基准风险率，在时点t时，某个影响因素就是使该基准风险率增至倍，即比例风险模型的基本结构为：
（11.9）
其中X为协变量向量，β为协变量的系数向量，，，是的转置向量。
2.模型拟合过程数据结构如下：
time
dead
treat
age
1
1
1
1
66
2
1
1
1
70
┆
┆
┆
┆
20
23
1
1
57
21
6
1
2
72
┆
┆
┆
┆
33
23
1
2
52
34
32
0
2
57
Stata命令：
use d:\data\cox
打开cox数据文件
gen ta=treat*age
生成age和treat的交互项ta
cox time treat age ta,dead(dead)
拟合含有交互项的Cox模型
cox命令的语句格式为：cox 生存时间变量 协变量,dead(结局变量)
输出结果如下：
①
Iteration 0,log likelihood = -68.781991
Iteration 1,log likelihood = -61.237091
Iteration 2,log likelihood = -59.671409
Iteration 3,log likelihood = -59.661862
Iteration 4,log likelihood = -59.661859
Refining estimates:
Iteration 0,log likelihood = -59.661859
Cox regression -- Breslow method for ties
Entry time 0 Number of obs = 34
② LR chi2(3) = 18.24
③ Prob > chi2 = 0.0004
Log likelihood = -59.661859 Pseudo R2 = 0.1326
------------------------------------------------------------------------------
time |
④
⑤
⑥
⑦
⑧
dead | Coef,Std,Err,z P>|z| [95% Conf,Interval]
-------------+----------------------------------------------------------------
treat | -3.005587 5.033301 -0.60 0.550 -12.87068 6.859502
age |,0891689,116237 0.77 0.443 -.1386514,3169892
ta |,0220285,0792198 0.28 0.781 -.1332394,1772965
------------------------------------------------------------------------------
输出结果中①为迭代过程；②为对模型H0假设：检验的结果，而本例中＝18.24得出相应P值为0.0004<<0.05，故拒绝H0，接受H1假设，认为不全为0；③为伪决定系数，描述了协变量的作用在总变异中所占的比例；④为对协变量系数的估计；⑤为协变量系数估计的标准误；⑥⑦为对协变量的检验，z近似服从标准正态分布，本例中先考察交互项ta，P>0.05，故认为treat和age之间没有交互作用；⑧为协变量系数的95%可信区间估计。
由于协变量之间没有交互作用，将交互项从模型中除去再进行分析。
cox time treat age,dead(dead)
拟合不含有交互项的Cox模型
输出结果为：
Iteration 0,log likelihood = -68.781991
Iteration 1,log likelihood = -59.857311
Iteration 2,log likelihood = -59.701191
Iteration 3,log likelihood = -59.700845
Refining estimates:
Iteration 0,log likelihood = -59.700845
Cox regression -- Breslow method for ties
Entry time 0 Number of obs = 34
LR chi2(2) = 18.16
Prob > chi2 = 0.0001
Log likelihood = -59.700845 Pseudo R2 = 0.1320
------------------------------------------------------------------------------
time |
dead | Coef,Std,Err,z P>|z| [95% Conf,Interval]
-------------+----------------------------------------------------------------
treat | -1.616596,4948552 -3.27 0.001 -2.586495 -.6466978
age |,119485,0407113 2.93 0.003,0396924,1992776
------------------------------------------------------------------------------
可见去除交互项后treat、ageP值均<0.05，在水平均有统计学意义；
至此，我们可以写出cox回归方程：
3.变量的危险比（risk ratio，记为RR）
RR=exp() （11.10）
表示协变量每增加一个单位，危险度改变多少倍。如协变量treat的= -1.617，RRtreat=0.199,表示treat变量水平2与1比较，treat=2的危险度是treat=1的0.199倍，提示治疗方案2优于治疗方案1。而age的=0.119，RRage=1.127，表明年龄每增加一岁，死亡的可能性增加到1.127倍。实际上，只需要在Stata命令中加入hr，就可以直接给出RR值。
cox time treat age,dead(dead) hr
输出结果：
------------------------------------------------------------------------------
time |
dead | Haz,Ratio Std,Err,z P>|z| [95% Conf,Interval]
-------------+----------------------------------------------------------------
treat |,1985735,0982651 -3.27 0.001,0752835,5237726
age | 1.126916,0458782 2.93 0.003 1.040491 1.220521
------------------------------------------------------------------------------
Haz.Ratio即为所需要的RR值。
至此，我们已经可以回答前面提出的3个问题了。
①在排除年龄的影响后，不同的治疗方案对生存时间的影响是不同的，方案2优于方案1；
②在排除不同治疗方案的影响后，年龄对生存时间的影响是不同的，年龄越大，越危险；
③没有理由认为治疗方案和年龄存在交互作用。
4.Cox模型分析注意事项：
（1）．默认假设为值随着自变量的变化而等比例变化，如本例中协变量age的值为0.119，即认为年龄每增加一岁，值增加0.119，而不论年龄是从55增加到56还是从70增加到71。
（2）．如果协变量不符合等比例假设的要求，那么就不能将连续性变量直接进入模型，而需要对其按段进行分层，必要时设置哑变量进行分析。
（3）．协变量既可以是连续变量，也可是是分类变量，还可以同时分析交互作用，这是与其他统计方法相比的重要优点。
（4）．选入模型的变量是统计学上的有关变量，有统计学意义并不一定都与生存时间有因果关系，需要结合具体案例进行分析。
思考题
1．明确生存率、生存概率、死亡率的关系。
2．两组淋巴瘤患者治疗后复发时间（月）如下，问两种治疗方案对淋巴瘤的缓解情况是否一样？
对照组
2
3
9
10
10
12+
15
15+
16
18+
24+
30
36+
40+
45+
处理组
9
12+
16
19
19+
20+
20+
24+
24+
30+
31+
34+
42+
44+
53+
3．28名自愿者参加了一项戒烟计划来帮助他们戒烟（原先所有的都是吸烟者），随访了35周，另外收集了性别、年龄、教育年数等因素，试分析这些因素和戒烟时间是否有关？
戒烟时间
week
是否戒烟失败
fail(1=fail)
性别
gender(1=male,2=female)
年龄
age
教育年数
edu
1
1
1
61
10
1
1
1
65
18
3
1
1
52
12
4
1
1
56
9
5
1
1
63
9
5
1
1
58
10
8
1
1
56
9
8
0
1
58
8
8
1
1
52
8
8
1
1
49
10
11
1
1
55
10
17
1
1
49
10
23
1
1
52
10
6
1
2
67
11
7
1
2
58
17
9
0
2
56
12
13
1
2
62
10
15
0
2
50
11
16
1
2
67
7
20
0
2
55
14
22
1
2
58
13
23
1
2
47
9
32
0
2
52
10
19
0
3
49
14
24
1
3
58
10
25
1
3
55
9
28
0
3
48
16
35
0
3
48
12
基本概念
1．1 生存时间、截尾例11.1：某医师从2002年1月1日起对６名肺癌患者进行跟踪观察，记录其结局，随访记录见表11.1：
表11.1 6例肺癌患者随访记录研究号
(1)
姓名
(2)
术后开始随访时间
(3)
终止随访时间
(4)
结局
(5)
生存天数
(6)
1
赵××
00-01-16
00-09-01
失去联系
229+
2
钱××
00-01-26
01-03-18
局部复发死亡
417
3
孙××
00-02-05
01-12-25
车祸死亡
689+
4
李××
00-02-19
02-10-01
研究终止
955+
5
王××
00-03-07
00-11-11
局部复发死亡
249
6
刘××
00-03-20
02-01-24
远处转移死亡
675
生存时间（survival time）是指从某起点事件(start point event)开始到被观测对象出现终点事件（endpoint event）所经历的时间。如从疾病“确诊”到“死亡”；从“治疗结束”到“复发”；某电子设备从“开始使用”到“出现故障”，等。由此可见，此处的“生存”是一个广义的概念。根据研究目的不同，可以有不同的“起点时间”和“终点时间”。生存时间的单位可以用年，月，周，日，甚至时，分，秒等表示。上表中的(6)即为生存时间，有2种类型：
（1）．完全数据（complete data）：是指被观测对象从观察起点至出现终点事件所经历的时间。在例11.1中，如果终点事件被设为死于肺癌，那么417、249、675为完全数据。
（2）．截尾数据（censored data）：是指在出现终点事件前，被观测对象的观测过程终止了。由于被观测对象所提供的信息是“不完全的”，只知道他们的生存时间超过了截尾时间。例11.1中，229+、689+、955+为截尾数据。
起点事件（start point event）：指研究者根据研究目的开始关心某一事件的起点，如上面所说的“确诊”、“治疗结束”、“开始使用”等等。
终点事件（endpoint event）：指研究者根据研究目的所关心的某一事件，如上面所说的“死亡”、“复发”、“出现故障”等等。需要注意一点，不同的研究目的有不同的终点事件，如果研究的是肿瘤的局部复发情况，那么死于肿瘤远处转移只能算做截尾，而不是终点事件。
截尾(censor)，又称终检、删失，主要有3种情况：①失访：指失去联系，如信访未回信，上门不见人，电话采访不答理或搬迁未留新地址等原因；②退出：指退出研究，如因其他非此次研究疾病死亡，临时改变治疗方案而中途退出研究；③终止：指研究时限已到而终止观察。
1．2 死亡概率、生存概率
死亡概率（mortality probability）：记为q，是指在某段时间开始时存活的个体在该段时间内死亡的可能性大小。若无截尾数据，死亡概率的估计公式为：
(11.1)
生存概率（survival probability）：记为p，是指在某段时间开始时存活的个体至该时间结束时仍存活的可能性大小。生存概率的公式为：
(11.2)
由于生存分析中常存在截尾值，假定失访等截尾事件在观察时段的各个时间点等机会发生的，故分母改用校正观察例数
1．3 生存率及标准误、生存曲线生存率（survival rate）：用S(tk)表示，是指病人经历tk个单位时间后仍存活的概率。若无截尾数据，则
(11.3)
其中t为病人的存活时间，但如果资料中含有截尾数据，分母必须按时段进行校正，此公式则不再适用，此时生存率的计算公式应为
(11.4)
其中、、…、表示不同时间段的生存概率，可以看出，生存率是多个时段生存概率的累积，故生存率又称为累积生存概率（cumulative probability of survival）。
生存率标准误的计算公式为
(11.5)
例11.2：某医院对100名恶性肿瘤术后病人进行随访所获得的资料如表11.2：
表11.2 某恶性肿瘤100例术后随访资料术后年数
t～
(1)
期初观察例数
n0
(2)
期内死亡例数
d
(3)
期内截尾例数
c
(4)
校正期初例数
nc=n0-c/2
(5)
死亡概率
q=d/n
(6)
生存概率
p=1-q
(7)
t+1年生存率
S(t+1)
(8)
生存率标准误
SE(S(t+1))
(9)
0～
100
10
0
100
0.1000
0.9000
0.9000
0.0300
1～
90
15
0
90
0.1667
0.8333
0.7500
0.0433
2～
75
16
8
71
0.2254
0.7746
0.5810
0.0501
3～
51
12
4
49
0.2449
0.7551
0.4387
0.0520
4～
35
14
7
31.5
0.4444
0.5556
0.2437
0.0484
5～
14
7
7
10.5
0.6667
0.3333
0.0812
0.0390
（1）为组段，（2）为各组段期初观察例数，（3）、（4）分别为各组段内的死于该恶性肿瘤例数和截尾例数，（5）为校正期初例数，（6）、（7）分别为各组段的死亡概率和生存概率，计算公式分别为(6)=(3)/(5)、(7)=1－(6)，（8）为t+1年生存率，（9）为t+1年生存率的标准误。
在本例中，q1=10/100，p1=1-q1=90/100，q2=15/90，p2=1-q2=75/90。
由式（11.4）求得2年生存率为S(2)= p1p2=(90/100)×(75/90) =75/100=0.75
由式（11.3）求得2年生存率为S(2)=(100－10－15)／100＝75/100=0.75
两式求得结果相同，但在继续计算S(3)、S(4)…时，由于式（11.3）不能对截尾数据进行处理，故不再适用，只能以（11.4）进行计算。
生存曲线（survival curve）：以时间为横轴，生存率为纵轴，将各个时点的生存率连接在一起的曲线图。
1．4 中位生存期中位生存期(Median Survival Time)：又称为半数生存期，即当累积生存率为0.5时所对应的生存时间，表示有且只有50%的个体可以活过这个时间。要注意中位生存期通常不等于生存时间的中位数（除非在这个时间点之前没有删失值存在）。
第二节 生存曲线估计
对于生存曲线的估计，我们介绍常用的2种方法：
1．Kaplan-meier法，又称乘积极限法（Product-Limit Method，简称PL法），由Kaplan-Meier在1958年提出，适用于样本量较小，难以将生存时间按组段划分，此时是利用tk时刻之前各时点上生存概率的乘积来估计在时刻tk的生存率，不需要对被估计的资料分布作任何假设。
例11.3：一组病人的生存时间（日）如下，用Kaplan-Meier法估计其生存曲线（+代表截尾）。
3,5,5,6+,8,16+,22,30,47+,71
为便于说明，将计算过程和结果列在表11.3，其步骤为：
（1）．将所有生存时间按从小到大排列（t），包括完全和截尾生存时间。
（2）．列出各期初暴露病例数（n），它是指在t时刻前仍存活的病例数。
（3）．将各期内死亡例数（d）和截尾例数（c）分别写在第（3）、（4）列。
（4）．计算各期的死亡概率q,q=d/n，截尾数不参加计算，结果见第（5）列。例如生存时间为5天时，q=2/9=0.2222。
（5）．计算各期的生存概率p,p=1-q，例如生存时间为5天时，p=1-0.2222=0.7778，计算结果见第（6）列。
（6）．计算各t时刻的生存率S(tk)。计算tk时刻生存率时可以用小于和等于tk时刻的各时点生存概率的乘积得到，计算结果见第（7）列。
表11.3 Kaplan-Meier法估计生存率计算表存活时间（天）
期初暴露病例数
期内死亡例数
期内截尾例数
死亡概率
生存概率
生存率
生存率标准误
t
n
d
c
q=d/n
p=1-q
S(tk)
SE(S(tk))
(1)
(2)
(3)
(4)
(5)
(6)
(7)
(8)
3
10
1
0
0.1000
0.9000
0.9000
0.0949
5
9
2
1
0.2222
0.7778
0.7000
0.1449
6
7
0
1
0.0000
1.0000
0.7000
0.1449
8
6
1
1
0.1667
0.8333
0.5833
0.1610
16
5
0
1
0.0000
1.0000
0.5833
0.1610
22
4
1
0
0.2500
0.7500
0.4375
0.1747
30
3
1
1
0.3333
0.6667
0.2917
0.1666
47
2
0
1
0.0000
1.0000
0.2917
0.1666
71
1
1
0
1.0000
0.0000
0.0000
.
用Stata软件计算时，数据库结构如下：
time
dead
1
3
1
2
5
1
3
5
1
4
6
0
5
8
1
6
16
0
7
22
1
8
30
1
9
47
0
10
71
1
Stata命令如下：
stset time,failure(dead=1)
将该数据库定义为生存资料数据库
sts list
输出各时点生存率
输出结果：
Beg,Net Survivor Std.
Time Total Fail Lost Function Error [95% Conf,Int.]
-------------------------------------------------------------------------------
3 10 1 0 0.9000 0.0949 0.4730 0.9853
5 9 2 0 0.7000 0.1449 0.3287 0.8919
6 7 0 1 0.7000 0.1449 0.3287 0.8919
8 6 1 0 0.5833 0.1610 0.2298 0.8207
16 5 0 1 0.5833 0.1610 0.2298 0.8207
22 4 1 0 0.4375 0.1747 0.1187 0.7256
30 3 1 0 0.2917 0.1666 0.0471 0.6085
47 2 0 1 0.2917 0.1666 0.0471 0.6085
71 1 1 0 0.0000,,,
-------------------------------------------------------------------------------
输出结果依次为时点、各时点的期初例数、死亡例数、截尾例数、生存函数（生存率）及其标准误、生存率的95%可信区间。如果需要输出生存函数曲线图，接着输入命令：
sts graph,title（Kaplan-Meier法对10名病人的生存率估计）
输出生存曲线
输出结果：
图11.1 Kaplan-Meier法对10名病人的生存率估计其中title()选项要求输出图像的标注。
2．寿命表法（Life-Table Method，简称LT法）是通过计数落入时间区间［tk-1,tk）内的失效和截尾的观察例数来估计该区间上的死亡概率，然后用该区间及其之前各区间上的生存概率之积来估计S(tk)。当样本量较大或者无法准确得知研究结果出现的时间时，可以将各研究对象的生存时间按某个时间段（年、月等）进行分组计算其生存率。
对例11.2用寿命表法计算生存率，具体计算步骤前面已有叙述，这里主要介绍stata的用法。
数据库结构如下：
time
num
dead
1
0
10
1
2
1
15
1
┆
┆
┆
┆
6
5
7
1
7
0
0
0
┆
┆
┆
┆
11
4
7
0
12
5
7
0
此时的数据库中有3个变量，time为生存时间（年），time=n表示该组患者在第n年初随访到了但是在n+1年初没有随访到。dead是观察结局，dead=1表示死亡，dead=0表示截尾。num为各组频数。
Stata程序如下：
use d:\data\ltable1
打开ltable1数据文件
ltable time dead [weight=num]
制作寿命表
输出结果：
Beg,Std.
Interval Total Deaths Lost Survival Error [95% Conf,Int.]
-------------------------------------------------------------------------------
0 1 100 10 0 0.9000 0.0300 0.8221 0.9449
1 2 90 15 0 0.7500 0.0433 0.6529 0.8236
2 3 75 16 8 0.5810 0.0501 0.4765 0.6718
3 4 51 12 4 0.4387 0.0520 0.3354 0.5371
4 5 35 14 7 0.2437 0.0484 0.1557 0.3425
5 6 14 7 7 0.0812 0.0390 0.0260 0.1779
-------------------------------------------------------------------------------
输出结果依次为各段生存时间起点与终点、期初人数、期内死亡人数、截尾值人数、生存率、标准误及其相应的95%可信区间。要输出生存函数曲线图，只须在ltable命令后加入graph选项：
ltable time dead [weight=num],graph notab noconf title (100例恶性肿瘤患者术后生存率情况)
图11.2 100例恶性肿瘤患者术后生存曲线本命令中，graph要求输出生存函数曲线；notab要求不输出寿命表；noconf要求在生存函数曲线中不显示95%可信区间，title（）则是给图形加上标注。
如果获得资料中每一例患者都有自己的生存时间和随访结局，并且样本量较大时，也可以用寿命表法进行计算，此时的stata数据库结构为：
time
dead
┆
┆
┆
stata命令为：
ltable time dead,interval(365)
以365天（1年）为组段制作寿命表
得到结果是一致的。
3．两种方法的比较：
①．寿命表法适用于大样本或无法准确得知研究结果出现时间的资料，Kaplan-Meier法主要用于小样本，也可以用于大样本。
②．寿命表法是按照指定的时段来分段，估计的是时间区间右端点上的生存率；Kaplan-Meier法是根据死亡时点分段，逐个估计死亡时点的生存率。
③．寿命表法不能确切得知死亡时间，假定每个时间段中的“死亡”是呈均匀分布，生存率为线性变化，故简单化以直线相连接；Kaplan-Meier法其生存曲线是左连续的阶梯型曲线，间断点的纵坐标在下一阶处，当样本量较大及死亡时点较多时，阶梯形就不明显了。
4．中位生存期的估计
寿命表法由于默认组段内生存率的变化是均匀的，因此可以直接在生存曲线上进行内插（图11.3a），如例11.2，3年生存率为0.5810，4年生存率为0.4387，当生存率为0.5000时，中位生存时间（年）为=3.57。
而在Kaplan-Meier法中，由于估计的是时点生存率，生存曲线是阶梯形的，对中位生存期的估计有两种观点：一种观点认为中位生存期为生存率降到0.5或以下的首个生存时间，另一观点认为需要先将生存率为0.5两侧左端点连线再进行内插。例11.3在两种观点下中位生存期分别为22天和16天（图11.3b）。stata软件默认的中位生存期估计为前一种观点，可通过stsum命令得到实现。
图11.3a 寿命表法中位生存期估计
图11.3b Kaplan-Meier法中位生存期估计
第三节 生存曲线的比较这一节将介绍两组生存资料比较最常用的时序检验（log-rank test），无效假设H0为两条总体生存曲线相同。如果H0成立，两组生存资料来自同一总体，用两组合并的资料估计一条生存函数曲线，该生存函数在各时段中所计算的理论死亡人数与实际死亡人数相差不会太大，否则拒绝无效假设，接受备择假设，认为各组总体生存曲线不同或不全相同。
例11.4：某医生收集23例晚期肺癌患者在接受化疗后的生存时间t（月），按接受治疗方案的不同划分为2组（1为常规方案，2为新方案），问不同的治疗方案对其生存时间长短的影响有无显著性差异。
常规方案组：
3
3+
4
6
6+
7+
8
9+
14
19
新方案组：
5
9
9+
10+
11
11
12
13+
17
18+
19+
20
32+
计算步骤：
（1）建立假设
H0：两总体生存函数曲线相同
H1：两总体生存函数曲线不同
（2）按时间排序将两组未截尾的完全生存时间从小到大混合排序，见表11.4第（2）列表11.4 两组肺癌患者资料分析表序号
生存时间
死亡数
期初观察数
理论死亡数
实际数－死亡数
j
tj
d1j
d2j
n1j
n2j
e1j
e2j
d1j-e1j
d2j-e2j
(1)
(2)
(3)
(4)
(5)
(6)
(7)
(8)
(9)
(10)
1
3
1
0
10
13
0.435
0.565
0.565
-0.565
2
4
1
0
8
13
0.381
0.619
0.619
-0.619
3
5
0
1
7
13
0.350
0.650
-0.350
0.350
4
6
1
0
7
12
0.368
0.632
0.632
-0.632
5
8
1
0
4
12
0.250
0.750
0.750
-0.750
6
9
0
1
3
12
0.200
0.800
-0.200
0.200
7
11
0
2
2
9
0.364
1.636
-0.364
0.364
8
12
0
1
2
7
0.222
0.778
-0.222
0.222
9
14
1
0
2
5
0.286
0.714
0.714
-0.714
10
17
0
1
1
5
0.167
0.833
-0.167
0.167
11
19
1
0
1
3
0.250
0.750
0.750
-0.750
12
20
0
1
0
2
0.000
1.000
0.000
0.000
合计
6
7
3.272
9.728
2.728
-2.728
（3）将不同生存时间的死亡数按组别归入（3）、（4），期初观察数按组别归入（5）、（6），计算时，先将两组每一生存时间tj的资料列成一个2×2表的形式，见表(11.5)，
表11.5 2×2表形式
死亡
生存
合 计
常规方案组
d1j
n1j-d1j
n1j
新方案组
d2j
n2j-d2j
n2j
合 计
Dj
Sj
Nj
然后按照四格表χ2 检验的方法计算理论死亡数。即
eij = ,i=1,2 (11.6)
将不同生存时间的理论死亡数计算结果列在（7）（8）。常规方案组的理论死亡数合计是3.272，新方案组的理论死亡数合计是9.728。全部实际死亡数与理论死亡数之差为：
O1 - E1 = 2.728
O2 - E2 =-2.728
由结果可看出，两组的实际数与理论数之差值是一样的，只是符号相反。因此计算时，选择任一组计算的结果就可以。
对于两组生存率进行Log Rank 检验，其统计量服从χ2 分布，计算公式为：
,＝1 (11.7)
其中
= ,i =1,2 (11.8)
本例中，=2.260，故
χ2 === 3.29
式（11.7）的近似公式为
（11.9）
本例中，3.04，结果略有不同。
下面用Stata软件进行计算：
数据库结构：
time
treat
dead
1
3
常规方案组
1
2
3
常规方案组
0
┆
┆
┆
┆
9
9
常规方案组
0
10
9
新方案组
1
┆
┆
┆
┆
22
20
新方案组
1
23
32
新方案组
0
Stata命令：
use d:\data\km2
打开km2数据库
stset time,failure(dead)
定义为生存数据库
sts list,by(treat)
按组别输出不同时点生存率
sts test treat,logrank
对两总体生存率曲线进行logrank检验
sts graph,by(treat)
按组别输出生存曲线
输出结果：
failure _d,dead
analysis time _t,time
Beg,Net Survivor Std.
Time Total Fail Lost Function Error [95% Conf,Int.]
-----------------------------------------------------------------------------
常规方案组
3 10 1 1 0.9000 0.0949 0.4730 0.9853
4 8 1 0 0.7875 0.1340 0.3809 0.9426
6 7 1 1 0.6750 0.1551 0.2906 0.8825
7 5 0 1 0.6750 0.1551 0.2906 0.8825
8 4 1 0 0.5063 0.1868 0.1396 0.7903
9 3 0 1 0.5063 0.1868 0.1396 0.7903
14 2 1 0 0.2531 0.2019 0.0138 0.6438
19 1 1 0 0.0000,,,
新方案组
5 13 1 0 0.9231 0.0739 0.5664 0.9888
9 12 1 1 0.8462 0.1001 0.5122 0.9591
10 10 0 1 0.8462 0.1001 0.5122 0.9591
11 9 2 0 0.6581 0.1407 0.3200 0.8576
12 7 1 0 0.5641 0.1488 0.2436 0.7928
13 6 0 1 0.5641 0.1488 0.2436 0.7928
17 5 1 0 0.4513 0.1560 0.1549 0.7121
18 4 0 1 0.4513 0.1560 0.1549 0.7121
19 3 0 1 0.4513 0.1560 0.1549 0.7121
20 2 1 0 0.2256 0.1776 0.0151 0.5896
32 1 0 1 0.2256 0.1776 0.0151 0.5896
-----------------------------------------------------------------------------
failure _d,dead
analysis time _t,time
Log-rank test for equality of survivor functions
| Events Events
treat | observed expected
------------+-------------------------
常规方案组 | 6 3.27
新方案组 | 7 9.73
------------+-------------------------
Total | 13 13.00
chi2(1) = 3.29
Pr>chi2 = 0.0696
图11.4 两组肺癌患者的生存曲线比较
从输出结果可以看出，logrank检验结果为=3.29，=1，P=0.0696，按水平，没有理由拒绝H0假设，认为两组的总体生存曲线相同，在本例中，即在晚期肺癌中，没有理由认为新化疗方案和常规化疗方案的疗效不同。
需要注意的是，在做Log-Rank检验的时候，除了生存资料的基本要求之外，还要求各组生存曲线不能交叉。因为这种交叉提示可能有混杂因素，应该采用分段分析或者采用多因素方法来分析。
第四节 指数模型和Weibull模型设生存时间T，如果生存率S(t)满足下列模型：
其中回归系数，协变量，
因此
称生存时间T服从指数分布模型，显然越大，也越大，因此生存率就越小。故称为指数回归模型的风险函数。称为指数回归模型的累积风险函数。
如果生存率S(t)满足下列模型：
则称生存时间T服从Weibull回归模型，显然当(＝1时，Weibull回归模型就是指数回归模型，因此指数回归模型是Weibull回归模型的特例。Weibull回归模型的累积风险函数为，相应的风险函数为。
以上可以归纳为累积风险函数为－ln(S(t))，并且Weibull模型的累积风险函数为－ln(S(t))＝表示累积风险以t的幂函数增长。特别指数分布模型的累积风险以t直线增长。
以下举一个例子说明这类模型的应用。
例11.5：某医生在研究2种治疗方案对肝癌术后患者的疗效，收集了60名患者的生存资料，随机分为两组，分别采用2种不同的治疗方案，同时收集了患者的年龄情况，生存时间的单位为年，死亡的结局定义为1，生存的结局定义为0，资料如下：
方案A(drug=0)
方案B(drug=1)
年龄
结局
生存时间
年龄
结局
生存时间
71
1
0.5
54
0
26.6
70
1
0.1
69
1
0.1
56
1
2.1
66
1
0.6
52
1
0.2
68
1
1.6
67
1
0.4
70
1
0.7
60
1
1.1
53
1
3.1
74
1
0.1
68
1
3.1
53
1
1.3
51
1
2.4
67
1
0.1
66
1
1
53
1
1.8
76
1
1.7
79
1
0.1
57
1
1.1
56
1
0.6
51
1
0.8
57
1
3.6
62
1
7.2
54
1
0.3
76
1
2.3
72
1
0.3
65
1
4.4
61
0
3.5
75
1
2.9
61
1
0.8
56
1
2.2
79
1
0.4
66
1
1.1
67
1
0.2
57
1
1.2
73
1
0.7
78
1
0.1
52
1
0.3
58
1
3.4
74
1
1.2
53
0
2.9
57
1
0.8
62
1
2.4
63
1
0.3
71
1
4
53
1
0.6
76
1
0.8
56
1
2.3
63
1
0.2
51
1
2
62
1
1.5
55
1
2
76
1
0.3
58
1
0.4
51
1
1.9
68
1
0.1
70
1
0.6
Stata数据格式如下
drug
age
d
t
0
71
1
0.5
0
70
1
0.1
0
56
1
2.1
0
52
1
0.2
0
67
1
0.4
0
60
1
1.1
0
74
1
0.1
0
53
1
1.3
0
67
1
0.1
0
53
1
1.8
0
79
1
0.1
0
56
1
0.6
0
57
1
3.6
0
54
1
0.3
0
72
1
0.3
0
61
0
3.5
0
61
1
0.8
0
79
1
0.4
0
67
1
0.2
0
73
1
0.7
0
52
1
0.3
0
74
1
1.2
0
57
1
0.8
0
63
1
0.3
0
53
1
0.6
0
56
1
2.3
0
51
1
2
0
55
1
2
0
58
1
0.4
0
68
1
0.1
1
54
0
26.6
1
69
1
0.1
1
66
1
0.6
1
68
1
1.6
1
70
1
0.7
1
53
1
3.1
1
68
1
3.1
1
51
1
2.4
1
66
1
1
1
76
1
1.7
1
57
1
1.1
1
51
1
0.8
1
62
1
7.2
1
76
1
2.3
1
65
1
4.4
1
75
1
2.9
1
56
1
2.2
1
66
1
1.1
1
57
1
1.2
1
78
1
0.1
1
58
1
3.4
1
53
0
2.9
1
62
1
2.4
1
71
1
4
1
76
1
0.8
1
63
1
0.2
1
62
1
1.5
1
76
1
0.3
1
51
1
1.9
1
70
1
0.6
操作如下：
定义生存分析数据集
stset t,f(d)
指数回归分析命令
streg drug age,d(e) nohr
No,of subjects = 60 Number of obs = 60
No,of failures = 57
Time at risk = 110.4000007
LR chi2(2) = 33.64
Log likelihood = -91.361189 Prob > chi2 = 0.0000
------------------------------------------------------------------------------
_t | Coef,Std,Err,z P>|z| [95% Conf,Interval]
-------------+----------------------------------------------------------------
drug | -1.230369,2670223 -4.61 0.000 -1.753723 -.7070145
age |,0672908,0157282 4.28 0.000,0364642,0981175
_cons | -4.043394 1.001847 -4.04 0.000 -6.006979 -2.07981
------------------------------------------------------------------------------
A方案drug=1与B方案drug=0的风险函数比(hazard ratio)
，
同理年龄增加一岁的风险函数比
也可以利用Stata直接计算得到：streg drug age,d(e)
No,of subjects = 60 Number of obs = 60
No,of failures = 57
Time at risk = 110.4000007
LR chi2(2) = 33.64
Log likelihood = -91.361189 Prob > chi2 = 0.0000
------------------------------------------------------------------------------
_t | Haz,Ratio Std,Err,z P>|z| [95% Conf,Interval]
-------------+----------------------------------------------------------------
drug |,2921849,0780199 -4.61 0.000,1731282,4931142
age | 1.069606,0168229 4.28 0.000 1.037137 1.103092
------------------------------------------------------------------------------
由于指数回归模型是Weibull回归模型的特例，所以通常用Weibull回归模型要优于指数回归模型。Stata命令如下
streg drug age,d(w) nohr
相应的结果如下
Weibull regression -- log relative-hazard form
No,of subjects = 60 Number of obs = 60
No,of failures = 57
Time at risk = 110.4000007
LR chi2(2) = 26.04
Log likelihood = -91.355743 Prob > chi2 = 0.0000
------------------------------------------------------------------------------
_t | Coef,Std,Err,z P>|z| [95% Conf,Interval]
-------------+----------------------------------------------------------------
drug | -1.24424,2984532 -4.17 0.000 -1.829198 -.6592825
age |,0680308,0172641 3.94 0.000,0341938,1018678
_cons | -4.091576 1.104093 -3.71 0.000 -6.255558 -1.927593
-------------+----------------------------------------------------------------
/ln_p |,0105058,1003806 0.10 0.917 -.1862365,2072482
-------------+----------------------------------------------------------------
p | 1.010561,1014407,8300773 1.230288
1/p |,9895491,0993315,8128179 1.204707
------------------------------------------------------------------------------
本例资料用Weibull回归模型进行参数估计与指数模型的参数估计非常接近，说明本例资料服从指数分布的回归模型。
指数分布的回归模型与Weibull分布的回归模型均属于生存分析中的比例风险参数模型，即：协变量的参数表达式与生存时间呈比例，更一般的比例风险模型为
其中回归系数，协变量，
因此
当X的各个分量全为0时，，故称h0(t)为基准风险函数(Baseline Hazard function)。
生存分析的模型可以分为参数模型和半参数模型，指数分布的回归模型和Weibull分布的回归模型都属于参数模型，这些模型可以直接估计生存率，但是要求资料服从对应的分布。这对实际研究的要求很高，往往难以确认。对于研究生存问题的危险因素而言，可以考虑下列的半参数模型：Cox模型。这种模型对资料的要求大大降低了。
第五节 Cox回归
由于log-rank仅能分析一个因素，因此，在同时分析2个或2个以上因素对生存时间影响的时候就显得无能为力了，这时我们就需要通过Cox比例风险模型来解决这些问题。以例11.5加以说明：
例11.5：某医生在研究2种治疗方案对晚期肺癌的影响，收集了34名患者的生存资料，同时收集了患者的年龄因素，资料如下：
方案1治疗组(N=20)
方案2治疗组(N=14)
生存月数
年龄
生存月数
年龄
1
66
6
72
1
70
6+
70
2
64
7
63
3
57
9+
61
4
61
10+
54
4
72
11+
66
5
68
13
67
5
63
15+
55
8
61
16
72
8+
63
19+
55
8
57
20+
60
8
54
22
63
11
55
23
52
11
60
32+
57
12
54
12
67
15
56
17
54
22
62
23
57
希望回答以下问题：
①在同样的年龄情况下，不同的治疗方案对生存时间有无影响？
②在同一治疗方案的情况下，年龄对生存时间有无影响？
③治疗方案和年龄有无交互作用？即在不同的年龄段中，不同的治疗方案对生存期的影响是否不同？
1.模型的结构
在Cox比例风险模型中，假定有一个基准风险率，在时点t时，某个影响因素就是使该基准风险率增至倍，即比例风险模型的基本结构为：
（11.9）
其中X为协变量向量，β为协变量的系数向量，，，是的转置向量。
2.模型拟合过程数据结构如下：
time
dead
treat
age
1
1
1
1
66
2
1
1
1
70
┆
┆
┆
┆
20
23
1
1
57
21
6
1
2
72
┆
┆
┆
┆
33
23
1
2
52
34
32
0
2
57
Stata命令：
use d:\data\cox
打开cox数据文件
gen ta=treat*age
生成age和treat的交互项ta
cox time treat age ta,dead(dead)
拟合含有交互项的Cox模型
cox命令的语句格式为：cox 生存时间变量 协变量,dead(结局变量)
输出结果如下：
①
Iteration 0,log likelihood = -68.781991
Iteration 1,log likelihood = -61.237091
Iteration 2,log likelihood = -59.671409
Iteration 3,log likelihood = -59.661862
Iteration 4,log likelihood = -59.661859
Refining estimates:
Iteration 0,log likelihood = -59.661859
Cox regression -- Breslow method for ties
Entry time 0 Number of obs = 34
② LR chi2(3) = 18.24
③ Prob > chi2 = 0.0004
Log likelihood = -59.661859 Pseudo R2 = 0.1326
------------------------------------------------------------------------------
time |
④
⑤
⑥
⑦
⑧
dead | Coef,Std,Err,z P>|z| [95% Conf,Interval]
-------------+----------------------------------------------------------------
treat | -3.005587 5.033301 -0.60 0.550 -12.87068 6.859502
age |,0891689,116237 0.77 0.443 -.1386514,3169892
ta |,0220285,0792198 0.28 0.781 -.1332394,1772965
------------------------------------------------------------------------------
输出结果中①为迭代过程；②为对模型H0假设：检验的结果，而本例中＝18.24得出相应P值为0.0004<<0.05，故拒绝H0，接受H1假设，认为不全为0；③为伪决定系数，描述了协变量的作用在总变异中所占的比例；④为对协变量系数的估计；⑤为协变量系数估计的标准误；⑥⑦为对协变量的检验，z近似服从标准正态分布，本例中先考察交互项ta，P>0.05，故认为treat和age之间没有交互作用；⑧为协变量系数的95%可信区间估计。
由于协变量之间没有交互作用，将交互项从模型中除去再进行分析。
cox time treat age,dead(dead)
拟合不含有交互项的Cox模型
输出结果为：
Iteration 0,log likelihood = -68.781991
Iteration 1,log likelihood = -59.857311
Iteration 2,log likelihood = -59.701191
Iteration 3,log likelihood = -59.700845
Refining estimates:
Iteration 0,log likelihood = -59.700845
Cox regression -- Breslow method for ties
Entry time 0 Number of obs = 34
LR chi2(2) = 18.16
Prob > chi2 = 0.0001
Log likelihood = -59.700845 Pseudo R2 = 0.1320
------------------------------------------------------------------------------
time |
dead | Coef,Std,Err,z P>|z| [95% Conf,Interval]
-------------+----------------------------------------------------------------
treat | -1.616596,4948552 -3.27 0.001 -2.586495 -.6466978
age |,119485,0407113 2.93 0.003,0396924,1992776
------------------------------------------------------------------------------
可见去除交互项后treat、ageP值均<0.05，在水平均有统计学意义；
至此，我们可以写出cox回归方程：
3.变量的危险比（risk ratio，记为RR）
RR=exp() （11.10）
表示协变量每增加一个单位，危险度改变多少倍。如协变量treat的= -1.617，RRtreat=0.199,表示treat变量水平2与1比较，treat=2的危险度是treat=1的0.199倍，提示治疗方案2优于治疗方案1。而age的=0.119，RRage=1.127，表明年龄每增加一岁，死亡的可能性增加到1.127倍。实际上，只需要在Stata命令中加入hr，就可以直接给出RR值。
cox time treat age,dead(dead) hr
输出结果：
------------------------------------------------------------------------------
time |
dead | Haz,Ratio Std,Err,z P>|z| [95% Conf,Interval]
-------------+----------------------------------------------------------------
treat |,1985735,0982651 -3.27 0.001,0752835,5237726
age | 1.126916,0458782 2.93 0.003 1.040491 1.220521
------------------------------------------------------------------------------
Haz.Ratio即为所需要的RR值。
至此，我们已经可以回答前面提出的3个问题了。
①在排除年龄的影响后，不同的治疗方案对生存时间的影响是不同的，方案2优于方案1；
②在排除不同治疗方案的影响后，年龄对生存时间的影响是不同的，年龄越大，越危险；
③没有理由认为治疗方案和年龄存在交互作用。
4.Cox模型分析注意事项：
（1）．默认假设为值随着自变量的变化而等比例变化，如本例中协变量age的值为0.119，即认为年龄每增加一岁，值增加0.119，而不论年龄是从55增加到56还是从70增加到71。
（2）．如果协变量不符合等比例假设的要求，那么就不能将连续性变量直接进入模型，而需要对其按段进行分层，必要时设置哑变量进行分析。
（3）．协变量既可以是连续变量，也可是是分类变量，还可以同时分析交互作用，这是与其他统计方法相比的重要优点。
（4）．选入模型的变量是统计学上的有关变量，有统计学意义并不一定都与生存时间有因果关系，需要结合具体案例进行分析。
思考题
1．明确生存率、生存概率、死亡率的关系。
2．两组淋巴瘤患者治疗后复发时间（月）如下，问两种治疗方案对淋巴瘤的缓解情况是否一样？
对照组
2
3
9
10
10
12+
15
15+
16
18+
24+
30
36+
40+
45+
处理组
9
12+
16
19
19+
20+
20+
24+
24+
30+
31+
34+
42+
44+
53+
3．28名自愿者参加了一项戒烟计划来帮助他们戒烟（原先所有的都是吸烟者），随访了35周，另外收集了性别、年龄、教育年数等因素，试分析这些因素和戒烟时间是否有关？
戒烟时间
week
是否戒烟失败
fail(1=fail)
性别
gender(1=male,2=female)
年龄
age
教育年数
edu
1
1
1
61
10
1
1
1
65
18
3
1
1
52
12
4
1
1
56
9
5
1
1
63
9
5
1
1
58
10
8
1
1
56
9
8
0
1
58
8
8
1
1
52
8
8
1
1
49
10
11
1
1
55
10
17
1
1
49
10
23
1
1
52
10
6
1
2
67
11
7
1
2
58
17
9
0
2
56
12
13
1
2
62
10
15
0
2
50
11
16
1
2
67
7
20
0
2
55
14
22
1
2
58
13
23
1
2
47
9
32
0
2
52
10
19
0
3
49
14
24
1
3
58
10
25
1
3
55
9
28
0
3
48
16
35
0
3
48
12
