11生存分析与SAS程序.docx
《11生存分析与SAS程序.docx》由会员分享，可在线阅读，更多相关《11生存分析与SAS程序.docx（31页珍藏版）》请在冰豆网上搜索。
11生存分析与SAS程序
第十一章生存分析与SAS程序
在医学研究中，考察研究因素对研究对象的效应，经典的研究设计方法只观察其所出现的结局或结果。
但是在一些诸如癌症等险恶疾病和慢性病的随访研究中，只观察其结局就不够全面，还需要考察研究对象出现某种结局所经历的时间，这段时间称为生存时间（survivaltime）。
生存分析是用来分析生存时间资料的统计方法，是近几十年来产生且发展甚为迅速的一门应用统计的分支。
到目前为止，生存分析已形成一套完整的体系，包括描述生存规律和进行组间比较的参数和非参数方法以及分析影响生存期因素的回归模型等。
本章主要介绍描述生存规律的寿命表和乘积限估计法，用于组间比较的logrank检验和作预后因素分析的Cox回归模型以及相应的SAS过程。
11.1生存分析的基本概念
11.1.1生存时间
在医学研究中，对于肿瘤、心血管等慢性疾病，要考察其治疗方法优劣，疾病预后的好坏以及影响疾病预后的因素，通常采用随访研究的方法。
对某一疾病作随访研究时，一般是从某一时间开始，观察到某一规定时间截止，而研究对象是始点以后陆续进入观察。
随访中要规定一个事件作为随访结局，例如：
病人死于研究疾病。
如果病人的随访结果是规定的结局，则称为失效或死亡（Failure），那么病人从进入观察随访到规定的结局出现，其间所经历的这段时间称为生存时间或失效时间（Survivaltime,Failuretime,Waitingtime）。
通常用Ti表示（i=1,2,…,n,n为观察个体数）。
在随访中，由于客观条件限制，不能也不可能将全部观察对象都观察到规定的结局（Failure）。
在观察截止时，一组观察对象除了出现规定的结局外，还有三种结果：
（1）死于其他疾病；
（2）由于迁移等原因失去联系（失访）；（3）随访截止时尚未出现规定的结局。
这三种结果虽原因不一，但提供的信息是一致的，这类病人不能获得确切的生存时间，而只知道其生存期比随访观察到的时间长的信息。
我们将此类病例数据称为截尾（Censored）数据。
其观察到的时间用Ci表示。
在获得生存时间时要注意以下几点：
1.随访结局可以不是死亡。
随访的结局根据研究问题的性质和目的的不同可以是死亡、复发、恢复等。
但要明确时间界线。
例如：
肿瘤病人已作了手术切除治疗，我们随访观察他有否复发，复发就是我们规定的结局。
2.确定进入随访观察的起点时间。
随访研究中，无论是同时进入观察还是陆续进入观察，对所有对象都要规定某一事件作为进入观察起点事件。
例如：
病人确诊，治疗开始或手术时间，动物试验中染毒等，作为起点事件要有明确的时间界线。
3.时间尺度。
医学研究一般用日历时间作为时间尺度。
但在某些医学问题研究中，可能不以日历时间作为时间尺度更能反映所研究的问题。
例如：
对儿童龋齿出现的随访研究中常以儿童出到第几颗牙才出现龋齿作为时间尺度。
11.1.2生存数据的特征
生存数据在结构上有它自己的特点，主要表现为：
1.生存时间数据的分布与常见的统计数据分布有明显不同，常呈指数分布、Weibull分布、对数正态分布、对数Logistic分布、Gamma分布或更为复杂的分布，生物医学中的生存时间分布有时呈现不规则状态，因而难以用传统的统计方法对这类数据进行处理。
2.生存数据中往往包含有截尾数据。
在随访研究中，若能观察到研究所规定的结局就能获得确切的生存时间，这类数据称为完全数据。
但在实际工作中，由于时间和客观条件限制，部分患者难以观察到终点事件，以致不能获得确切的生存时间，而只能获得进入观察至失访时这段时间。
在分析总结时若剔除这部分对象往往导致样本有偏性，并且这部分数据仍然提供了其实际生存时间大于观察到的截尾时间的信息，应充分利用这部分资料的信息。
因此，一般随访研究的样本数据中不可避免地包含有截尾数据。
11.1.3生存时间函数
描述生存时间分布规律的函数主要有生存函数、死亡概率函数、概率密度函数和危险率函数。
为了后文叙述方便，这里主要介绍生存函数和危险率函数。
1.生存函数
在描述生存规律的数量指标中，以往常用的指标是某个特定时间的生存率（例如：
3年生存率、5年生存率）。
这一指标的主要缺陷为不能反映整个生存规律，一个理想的指标应该是任意时间的生存率，即生存率是任意时刻t的函数。
其意义是研究个体生存时间长于t的概率。
若令T为生存期，s（t）为任意时刻t的生存率，得
s（t）=p（T≥t）0≤t<∞（11.1）
则称s（t）为生存率函数，简称生存函数。
（SurvivorFunction）。
将s（t）对t在直角坐标系作图画出生存率曲线，可从图上粗略估计出中位生存期或半数生存期。
即生存曲线上取生存率为50%时所对应的时间。
2.危险率函数（HazardsFunction）
如果我们考虑已活到t时刻的患者，在时间t附近的瞬间死亡危险性，根据数学上极限性质，可表示为：
h（t）=
（11.2）
则称h（t）为危险率函数。
即相当于条件瞬间死亡率或年龄别死亡率。
这是生存分析中一个特有的函数，为处理生存分析问题带来了极大方便。
3.s（t）与h（t）的关系
s（t）=
=
，其中H（t）=
，称H（t）为累计危险率函数。
生存函数s（t）和危险率函数h（t）在不同的生存时间分布中（例如Weibull、loglogistic等）有着特定的函数形式。
这里不再赘述，请参考有关专著和文献。
11.1.4生存分析的主要任务
1.根据一组特定的生存数据，对生存率函数s（t）、危险率函数h（t）以及中位生存期等指标进行估计。
常用的估计方法为寿命表估计法和乘积限估计法（product-limitestimates，又称Kaplan-Meier法）。
2.对不同的观察组，检验s（t）有无统计学意义，从而（举例说）比较不同治疗方法或药物的临床效果。
常用的单因素组间比较的方法有对数秩和检验（logranktest），广义Wilcoxon检验。
3.对可能影响生存期或预后的各因素在各种回归模型下进行变量筛选和评价。
常见的生存分析模型有Cox回归模型。
4.为达到以上目的，进行恰当的流行病学和临床上的设计，例如：
确定样本含量的大小、Failure数据的比例和观察截止时间等。
11.2生存率估计与组间比较
11.2.1生存率估计
生存率估计方法有参数法和非参数法两类。
参数法是根据专业知识和实际问题本身的特性或应用一定的统计手段选择特定的生存时间分布形式（如前面提到的指数分布、Weibull分布、对数正态分布等）拟合实际资料，求得特定分布中的参数，以此来描述生存和死亡规律。
这种方法计算较为复杂，但若资料确实服从某种特定的分布，那末对资料内在的特点和规律的表达更为紧凑和精炼。
非参数法对生存资料的分布型没有相应的要求，适用面比较广。
医学研究中，大量的生存资料其分布是不规则、不确定或未知分布，因而常用非参数法估计生存率。
根据样本含量的大小可分别选择寿命表法或乘积限估计法。
1.寿命表
寿命表法适用于样本例数较多的生存分析资料。
其计算生存率的基本原理和过程是首先将整个随访时间划分成若干个时间区间，然后分别计算每个时间区间开始时的观察病人数、死亡数和失访数，进而计算每个时间区间的死亡概率和生存概率。
随访开始至t时刻时尚生存也就相当于随访对象在t时间前各个时间区间均生存。
那末，根据概率的乘法法则，t时刻的生存率为t时刻前各时间区间生存概率的乘积。
下面结合实例简要介绍寿命表计算生存率的步骤。
例11.1某医院随访观察296例晚期肝癌患者确诊后的生存情况，得表11.1，试计算其确诊后的逐月生存率。
（1）确定分组区间ti-。
根据随访时间和病例生存时间的长短以及随访病例的多少确定分组的多少和区间宽度，最后一个区间的终点是在无穷大处。
本例以月为区间宽度，分成12个时间区间。
（2）计算期初观察人数ni，期内死亡数di和失访数wi。
首先计算第一个时间区间t0-的n0、d0和w0，然后计算下一个区间的n1，n1=n0-d0-w0，以此类推，得：
ni+1=ni-di-wi（11.3）
表11.1296例晚期肝癌患者生存率计算（时间单位：
月）
随访月数
（1）
期初观察人数
（2）
期内死亡人数
（3）
期内失访人数
（4）
校正人数
（5）
死亡概率
（6）
生存概率
（7）
生存率
（8）
生存率的标准误
（9）
0-
296
94
10
291
0.3230
0.6770
0.6770
0.0274
1-
192
74
15
184.5
0.4011
0.5989
0.4055
0.0294
2-
103
22
10
98
0.2245
0.7755
0.3144
0.0285
3-
71
22
6
68
0.3235
0.6765
0.2127
0.0263
4-
43
5
5
40.5
0.1235
0.8765
0.1864
0.0255
5-
33
6
6
30
0.2000
0.8000
0.1492
0.0246
6-
21
4
1
20.5
0.1951
0.8049
0.1201
0.0237
7-
16
2
1
15.5
0.1290
0.8710
0.1046
0.0230
8-
13
3
2
12
0.2500
0.7500
0.0785
0.0217
9-
8
2
0
8
0.2500
0.7500
0.0589
0.0202
10-
6
2
2
5
0.4000
0.6000
0.0353
0.0177
11及以上
2
2
—
2
1.0000
0.0000
0.0000
—
（3）校正观察人数ni’的计算由于在ti-区间内的失访人数wi在该区间内并未观察至区间的终点，因而该区间内的有效人数并非为ni。
设wi个失访个体在该区间内为均匀分布，即假定wi个体平均观察半个区间，那末校正人数为:
ni’=ni-wi/2（11.4）
（4）死亡概率qi和生存概率pi计算qi是指在ti-区间始端仍活着的病人在该区间内死亡的可能性，其计算公式为:
qi=di/ni′（11.5）
本例qi计算见表11.1第（6）栏。
pi是指在ti-区间始端仍活着的病人在该区间内生存的概率，其计算公式为:
pi=1-qi（11.6）
本例计算结果见表11.1第（7）栏。
（5）生存率s（ti）计算。
根据概率的乘法法则，s（ti）的计算公式为:
s（ti）=p0·p1·p2…pi-1（11.7）
本例计算结果见表11.1第（8）栏
（6）生存率标准误SE[s（ti）]的计算。
生存率是依据实际资料（样本）计算得到的，故有抽样波动，须对其标准误SE[s（ti）]进行估计，其近似计算公式为:
SE[s（ti）]=s（ti）·
（11.8）
求得生存率和标准误，就可以对这批随访对象在整个随访期间的生存规律进行描述；并可按正态近似原理，对总体生存率的可信区间作出估计。
2.乘积限估计法
当数据个体较少时，为充分利用每个数据的信息，必须采用更为精确的估计方法。
这些估计方法中应用最多、效率较高的是Kaplan-Meier在1958年提出的乘积限估计（Product-limitestimator）。
因而此法又称Kaplan-Meier法。
该法的基本思想与寿命表法基本相同，所不同的是将生存时间（包括截尾数据）逐个由小到大依次排列，并对其中的每个死亡点进行死亡概率、生存概率和生存率进行估计。
下面结合实例对其计算方法作简要介绍。
例11.2某中医研究院用猪苓提取物治疗患白血病患者16例,其生存数据经有序排列后如下,计算其不同时间的生存率。
2+，4，6+，6+，7.5+,8.5,9+,10,12+,13,18,19+,24,26,31,43+
（1）对生存时间进行编号和排序。
将生存时间（包括截尾时间）从小到大排序并编号。
遇相同生存时间只排一个；但有相同截尾时间或生存时间与截尾时间相同时，则分别列出。
本例见表11.2
（1）、
（2）栏。
表11.216例急性白血病患者生存数据比乘积限估计
顺号i
（1）
生存时间
（2）
观察人数
（3）
死亡人数
（4）
死亡概率
（5）
生存概率
（6）
生存率s（ti）（7）
标准误
（8）
1
2
16
0
0
1
1
2
4
15
1
0.0667
0.9333
0.9333
3
6
14
0
0
1
0.9333
0.0644
4
6
13
0
0
1
0.9333
5
7.5
12
0
0
1
0.9333
6
8.5
11
1
0.0909
0.9091
0.8485
0.0999
7
9
10
0
0
1
0.8485
8
10
9
1
0.1111
0.8889
0.7572
9
12
8
0
0
1
0.7542
0.1257
10
13
7
1
0.1429
0.8571
0.6464
0.1468
11
18
6
1
0.1667
0.8333
0.5387
0.1569
12
19
5
0
0
1
0.5387
13
24
4
1
0.2500
0.7500
0.4040
0.1657
14
26
3
1
0.3333
0.6667
0.2693
0.1558
15
31
2
1
0.5000
0.5000
0.1347
0.1231
16
43
1
0
0
1
（2）计算观察人数ni和死亡人数di。
先分别列出各死亡点上的死亡人数，然后计算每个死亡点或截尾点上的观察人数ni，其计算公式为:
ni+1=ni-di（wi）（11.9）
本例ni和di计算列于表11.2的（3）、（4）栏。
（3）计算死亡概率qi和生存概率pi。
qi和pi的计算公式为：
qi=di/ni（11.10）
pi=1-qi（11.11）
所有截尾点上的qi为0，pi为1。
本例qi和pi的计算列于表11.2的（5）与（6）栏。
（4）计算生存率s（ti）。
类似于寿命表，s（ti）的计算公式为：
s（ti）=p0·p1·p2…pi（11.12）
（5）生存率标准误SE[s（ti）]的计算，由于s（ti）是依据样本资料计算的，有抽样误差。
乘积限估计的标准误的近似计算公式为：
SE[s（ti）]=s（ti）·
（11.13）
本例s（ti）及标准误的计算列表11.2的（7）、（8）栏中。
11.2.2生存率的组间比较
在医学随访研究中，通常将病人按随机化方法分配到两种或多种治疗组中，然后随访观察和比较其生存时间的长短和生存率的大小，以此来考察各种治疗方案的优劣；或者分析和比较同一治疗方案下具有不同特征病人的生存率的大小，以此来探讨影响这种疗法的因素。
因此，生存率组间比较实际上是两条或多条生存曲线的比较。
生存率的假设检验方法有参数法和非参数法两类。
参数法要求生存时间已知服从于某种概率分布，对实际资料拟合分布并求得其相应的参数，然后通过比较不同组的分布参数来比较生存率是否相同。
非参数法对资料的分布没有要求，适用面比较广。
常见的有Logrank检验、Wilcoxon检验（Gehan检验）和似然比检验，似然比检验要求资料服从指数分布才有效。
这里主要介绍Logrank检验和Wilcoxon检验两种方法。
1.对数秩和检验（logranktest）
Logrank检验是Mantel等人在1966年提出的，这种方法是在组间生存率相同的检验假设（H0）下，对每组生存数据依据在各个时刻尚存活的患者数和实际死亡数计算期望死亡数，然后将期望死亡数与实际死亡数进行比较，作假设检验。
这种方法可适合两组或多组生存率比较。
这种方法在两组生存率比较时，计算比较简单，故结合实例对两组比较的计算方法作为详细介绍，然后推广到多组检验方法。
例11.3某中医研究院用猪苓提取物治疗急性的白血病患者16例（例11.2）,并设对照组10例，试比较两种疗法的优劣。
对照组的生存数据如下
对照组1.5,2+，3.5,6,6.5，6.5+，11,11+，13,17
（1）建立检验假设（H0）两种疗法治疗急性的白血病患者的生存率相同。
（2）将两组资料混合后统一排序并给定编号，当遇到生存时间与截尾时间相同时，只排一个时间点，但当相同数一个为生存时间，另一个为截尾时间，应分别列出且截尾时间在后。
（3）分别计算各ti时间上每一组的观察人数ngi，死亡数dgi和截尾wgi，以及两组合计的观察人数ni和死亡数di。
本例计算结果分别列于表11.3中的
（2）、（3）、（4）栏，（7）、（8）、（9）栏和（12）、（13）栏中。
（4）计算各组在时间ti上期望死亡数egi。
依据H0：
s1（t1）=s2（t2）,egi的计算公式为：
egi=ngi·di/ni（11.14）
（5）求各组实际和期望死亡数的和，并计算χ2统计量。
各组实际死亡数的和为Dg=
期望死亡数的和为Eg=
（l≤n）。
近似χ2统计量计算公式为：
χ2=
（11.15）
在H0条件下，χ2值服从自由度为1的χ2分布。
本例的χ2统计量为：
χ2=
+
=5.684
由于χ2=5.684>χ20。
05，1=3.84P<0.05，两组样本生存率间差异有统计学意义。
可以认为用猪苓提取物合用化疗治疗急性白血病的生存期大于对照疗法。
表11.3对数秩和检验计算表
编号
时间
用猪苓组（第1组）
对照组（第2组）
合计
i
ti
n1i
d1i
w1i
e1i
V1i
n2i
d2i
w2i
e2i
V2i
ni
di
（1）
（2）
（3）
（4）
（5）
（6）
（7）
（8）
（9）
（10）
（11）
（12）
（13）
1
1.5
16
0
0
0.615
0.2367
10
1
0
0.385
0.2367
26
1
2
2+
16
0
1
0.0
0.0
9
0
1
0.0
0.0
25
0
3
3.5
15
0
0
0.652
0.2268
8
1
0
0.348
0.2268
23
1
4
4
15
1
0
0.682
0.2169
7
0
0
0.318
0.2169
22
1
5
6
14
0
0
0.667
0.2381
7
1
0
0.333
0.2381
21
1
6
6+
14
0
2
0.0
0.0
6
0
0
0.0
0.0
20
0
7
6.5
12
0
0
0.667
0.2222
6
1
0
0.333
0.2222
18
1
8
6.5+
12
0
0
0.0
0.0
5
0
1
0.0
0.0
17
0
9
7.5+
12
0
1
0.0
0.0
4
0
0
0.0
0.0
16
0
10
8.5
11
1
0
0.733
0.1956
4
0
0
0.267
0.1956
15
1
11
9+
10
0
1
0.0
0.0
4
0
0
0.0
0.0
14
0
12
10
9
1
0
0.692
0.2130
4
0
0
0.308
0.2130
13
1
13
11
8
0
0
0.667
0.2222
4
1
0
0.333
0.2222
12
1
14
11+
8
0
0
0.0
0.0
3
0
1
0.0
0.0
11
0
15
12+
8
0
1
0.0
0.3025
2
0
0
0.0
0.0
10
0
16
13
7
1
0
1.556
0.1224
2
1
0
0.444
0.3025
9
2
17
17
6
0
0
0.857
0.0
1
1
0
0.143
0.1224
7
1
18
18
6
1
0
1.000
0.0
0
0
0
0.0
6
1
19
19+
5
0
1
0.0
0.0
0
0
0
0.0
5
0
20
24
4
1
0
1.000
0.0
0
0
0
0.0
4
1
21
26
3
1
0
1.000
0.0
0
0
0
0.0
3
1
22
31
2
1
0
1.000
0.0
0
0
0
0.0
2
1
23
43+
1
0
1
0.0
0
0
0
0.0
1
0
合计
8
11.788
2.1964
7
3.212
2.1964
精确法χ2统计量的计算公式为：
χ2=（Dg-Eg）2/vg（11.16）
式中vg为第g组的Fg的方差估计值，其计算公式为：
vg=
（11.17）
所计算的v1=v2，精确法χ2值服从自由度为1的χ2分布。
本例先求vgi，再求vg=
，得v1=v2=2.1964，χ2统计量为：
χ2=
=6.533专业结论同近似法。
对于多组生存率比较，其近似法χ2统计量的计算方法与两组相同。
设有m组g=1,2,…,m这时自由度υ=m-1。
关于多组Logrank检验的精确法χ2统计量计算较为复杂，其计算公式为：
χ2=sˊv-1sυ=m-1（11.18）
其中sˊ=（s1,s2,…,sg,…,sm-1）,sˊ为向量s的转置。
Sg的计算公式为：
sg=
=Og-Eg，g=1,…,m-1（11.19）
V为（m-1）（m-1）矩阵，记为V=[Vgh]（m-1）（m-1）。
其中Vgh为第g组与第k组之间的方差与协方差，计算公式为：
Vgh=
（11.20）
当g=h时，δgh=1；当g≠h时，δgh=0。
从多组生存率比较的计算公式可知，当m=2时,即为两组生存率的计算公式。
升级会员 