第1章 基本概念
第1节 生存资料的特点
生存资料(Survival Data)或失效时间资料(Failure-time Data)与多元线性回归资料很相似,只不过因变量(或反应变量)通常为观测对象生存的时间,常用t来表示。当然,生存时间是广义的,可以指在通常意义下生物体的生存时间、也可以指所关心的某现象(如疾病治愈后、合格品使用后)持续的时间。若生存时间是准确观测到的,则称为完全数据。生存资料的一个明显特点是:所收集的资料中常常包含不完全数据,也称为截尾数据、删失数据、终检数据(Censored Data)。包括删失数据的资料,称为删失资料。对于删失数据,既不能简单地弃之,踊能像对待完全数据那样给予充分的信任,需要采取一些技术处理。专门处理这种资料的统计方法,称为生存分析(Survival Analysis)。
导致数据删失有多种原因,最常见的有:失访(病人因搬家、随访信件丢失、车祸等原因,导致医生对他们的随访观察中断)和研究截止。由随机因素引起的,称为随机删失;若事先就定了截止日期,则称为定时删失(也称Ⅰ型删失);若事先就定了观察完多少例就
截止研究,则称为Ⅱ型删失(也称为定数删失)。在表达删失数据时,常在其右上角放一个“+”号;而用SAS软件分析时,常在其前放一个“-”号或产生1个指示变量(如:C=0表示删失数据、C=1表示完全数据,反过来也可以),便于计算时区别对待。为了使数据的表达与计算在形式上统一起来,本篇一律用负数表示删失数据,因生存时间不可能为负值,故不会产生混淆。
第2节 生存时间函数
描述生存时间规律的函数很多,统称为生存时间函数。其中最主要的有生存函数、死亡概率函数、概率密度函数和危险率函数。
1.生存函数(Survival Function)
生存函数也称为生存概率或累积生存率,常用S(t)表示,它表示一个体生存时间长于t的概率。在具体问题中,该函数在t时刻的取值可用式(5.1.1)来估计∶
S(t)≈生存时间长于t的病人数/病人总数 (5.1.1)
2.死亡概率函数(Failure Probability Function)
死亡概率函数简称为死亡概率,常用F(t)表示,它表示一个体从开始观察起到时刻t为止的死亡概率。它可以通过S(t)求得(详后)。
3.概率密度函数(Probability Density Function)
概率密度函数简称为密度函数,常用f(t)表示,它表示一个体死于(t,t+△t)小区间内的概率的极限。在具体问题中,该函数在t时刻的取值可用式(5.1.2)来估计∶
f(t)≈t时刻开始的区间内死汀人数/(病人总数×区间宽度) (5.1.2)
4.危险率函数(Hazard Function)
危险率函数也称为风险函数、瞬时死亡率、年龄别死亡率、条件死亡率,常用h(t)表示,它表示已存活到t的一个体,死于(t,t+△t)小区间内的概率的极限。在具体问题中,该函数在t时刻的取值可用式(5.1.3)来估计∶
h(t)≈t时刻开始的区间内死汀人数/(生存到t的病人数×区间宽度) (5.1.3)
5.上述几个函数之间的相互关系
(5.1.4)
(5.1.5)
(5.1.6)
(5.1.7)
上述各函数中“'”代表对t求导数,“∫”代表积分。
第3节 生存分析方法的分类
像普通统计分析一样,生存分析也有一套完整的方法:统计描述(包括求生存时间的分位数、中数生存期、平均数、生存函数的估计、判断生存时间的图示法);非参数检验(检验分组变量各水平所对应的生存曲线是否一致,常用的方法有对数秩检验(Log-rank Test)、威尔科克森检验(WilcoxonTest)和似然比检验(Likelihoodratio Test));COX模型(半参数模型)回归分析(在特定的假设之下,建立生存时间随多个危险因素变化的回归方程);参数模型回归分析(已知生存时间服从特定的参数模型时,拟合相应的参数模型,更准确地刻划变量之间的变化规律)。
第2章 生存资料的非参数统计方法
第1节 统计描述与非参数分析概述
1.统计描述
常用来反映一组生存时间平均水平的统计指标有中位数、平均数2种,因生存资料多为正偏态,故往往选用中位数更符合资料的特点。
对于寿命资料,首先需给出各时间点上生存函数的估计值,常用的方法有:乘积─极限法(Product-Limit Method,简称PL法)和寿命表法(Life-Table Method,简称LT法)。PL法是利用ti时刻之前各时间点上生存率的乘积来估计在时刻ti的生存函数S(ti)、而LT表法是通过计数落入时间区间[ti-1,ti]内的失效和删失的观察例数来估计S(ti)。
若能知道寿命函数的具体,可有的放矢地去选用相应的参数模型拟合资料,是非常有益的。实现这一目的途经是图解法,如:用(t, -logS(t))画图,若成一条直线,表明S(t)呈指数;又如:用(logt,log(-logS(t))画图,若成一条直线,表明S(t)呈图尔。当然,也有一些统
计检验方法,如:判断是否服从指数的G检验法、判断是否服从图尔的Mann-Scheuer-Fertig Tiku检验法和判断是否服从对数正态的W检验法等,具体检验方法参见有关专著。
2.各层间生存曲线的齐性检验
设全部受试者接受了k只同的处理,这k种处理实际上就是一个名义分类变量或枫因素的k个水平,于是,可按层估计生存函数。研究者常需比较k条生存曲线之间是否有显著差别,其方法有多种,SAS中用了以下3种:对数秩检验(Log-rank Test)、威尔科克森检验(Wilcoxon Test)和似然比检验(Likelihood ratio Test)。用它们来实现各层之间的齐性检验。
3.上述3种非参数检验的比较
当生存时间的为指数、图尔或属于比例危险模型时,Log-rank检验效率较高;当生存时间的为对数正态等时,Wilcoxon检验效率较高;似然比检验是建立在指数模型上的,故当资料偏离此模型时,其结果不如前2种检验方法稳健。
4.协变量与生存时间联系密切程度的检验
当资料中还包含与生存时间有关的其他连续变量(即协变量)时,也可分析它们与生存时间联系的密切程度。为实现此检验,LIFETEST过程中提供了2个分别建立在指数得分和
威尔科克森得分基础之上的删失数据线性秩统计量─Log-rank Test和Wilcoxon Test,这2种检验通过合并枫变量后进行计算,从而,校正了枫变量的影响。除了对重复(ties)生存时间的处理方法不同外,这里所说的2种检验与实现各层之间齐性检验中所提到的前2种检验是相同的。
为了不把读者的注意力引向复杂的计算,特将上述各种方法的具体计算公逝在本章第3节中再介绍,以便必要时备查。
第2节 用LIFETEST过程实现统计计算
[例5.2.1] 某医生收集到35例白血病患者治疗后的生存时间t(月),仔细观察后发现这些病人中有一部分人出现了白细胞(WBC)倍增的现象。现将他们按是否出现WBC倍增分成2组如下(注:负值代表删失数据),试用生存分析方法分析患者有无WBC倍增,对其生存时间长短有无显著影响。自艾自怨
A组(有WBC倍增):2,-2.5,3.5,4,4,-5,6,-6,7,-7,8,-9,10.5,12.5,19;
B组(无WBC倍增):2.5,5,7,-8.5,9,-10,11,-11,12,13,-14,15,-16,17,-18,19,-20,21,
24,32。
[SAS程序]──[D5P1.PRG]
DATA abc; PROC LIFETEST METHOD=PL
INFILE 'a:hlwbc.dat'; PLOTS=(S,LS,LLS);
INPUT lt @@; TIME t*censor(1);
IF lt<0 THEN censor=1; STRATA group;
ELSE censor=0; RUN;
IF _N_<16 THEN group='high-wbc'; PROC LIFETEST METHOD=LIFE
ELSE group='low-wbc'; PLOTS=(S,H);
t=ABS(lt); TIME t*censor(1);
STRATA group;
RUN;
(程序的第1部分) (程序的第2部分)
[程序修改指导] 用全部35个数据建立的数据文件名为HLWBC.DAT,第1个IF语句产生
1个指示变量CENSOR,其取值为1时为删失数据、取值为0时为完全数据。第2个IF语句产生1个分组变量GROUP,前15个数据属于有WBC倍增组、后20个数据属于无WBC倍增组。对表示删失和完全数据的变量lt取绝对值是为了保证参与计算的生存时间t都是正值。
第1个过程步是选择PL法计算(它也是隐含的方法)、第2个过程步是选择LT法计算。PLOTS= 要求绘图,其中S表示生存函数、L表示取对数、H表示危险率函数,图形的横坐标与纵坐标分别为:
S─(t,S)、LS─(t,-log(S))、LLS─(log(t),log(-(log(S))))、H─(t,H)
生存时间t与指示变量以乘法的形式写在TIME语句中、分组变量写在STRATA语句中。
当用寿命表(LT)法分析资料时,程序会自动形成生存时间的区间,也可人为指定生存时间的分组区间。做法是:在PROC语句的分号之前加上INTERVALS=(a to b by c),a、b、c分别为初值、终值、步长(必须是具体数值),步长的缺省值为1。
如果资料中还含有数值型的协变量,可将它们写在TEST语句中,如:TEST x1 x2 x3;以便检验协变量与生存时间联系的密切程度。当然,若有PHREG和LIFEREG过程,用它们建立起因变量t随自变量(即危险因素)变化的回归模型,可更好地揭示变量之间的内在
联系。
[输出结果及其解释]
Product-Limit Survival Estimates
GROUP = high-wbc
① ② ③ ④ ⑤ ⑥
Survival
Standard Number Number
T Survival Failure Error Failed Left
0.0000 1.0000 0 0 0 15
2.0000 0.9333 0.0667 0.0644 1 14
2.5000* . . . 1 13
3.5000 0.8615 0.1385 0.0911 2 12
4.0000 . . . 3 11
4.0000 0.7179 0.2821 0.1198 4 10
5.0000* . . . 4 9
6.0000 0.6382 0.3618 0.1304 5 8
6.0000* . . . 5 7
7.0000 0.5470 0.4530 0.1400 6 6
7.0000* . . . 6 5
8.0000 0.4376 0.5624 0.1487 7 4
9.0000* . . . 7 3
10.5000 0.2917 0.7083 0.1550 8 2
12.5000 0.1459 0.8541 0.1290 9 1
19.0000 0 1.0000 0 10 0
* Censored Obrvation
Quantiles 75% 12.5000 Mean 9.0775
50% 8.0000 Standard Error 1.6768
25% 4.0000
这是用PL法对第1组生存资料进行统计描述的结果。标号①~⑥分别是生存时间、生存
概率、死亡概率、生存概率的标准误差、已观察到的不同失效时间的数目、尚未观察到的不同失效或删失时间的数目,打*号的是删失观察值。接着,给出了生存时间的四分位数、均数及其标准误差。结果显示∶第1组患者中有25%的人(约4人)的生存时间短于4个月,即有75%的人的生存时间长于4个月;同理,可解释P50=8(个月)、P75=12.5(个月)的含义。由此可知∶该组患者的中数生存期为8个月、平均生存期约为9个月。
Product-Limit Survival Estimates
GROUP = low-wbc
Survival
Standard Number Number
T Survival Failure Error Failed Left
0.0000 1.0000 0 0 0 20
2.5000 0.9500 0.0500 0.0487 1 19
5.0000 0.9000 0.1000 0.0671 2 18
7.0000 0.8500 0.1500 0.0798 3 17
8.5000* . . . 3 16
9.0000 0.7969 0.2031 0.0908 4 15
10.0000* . . . 4 14
11.0000 0.7400 0.2600 0.1006 5 13
11.0000* . . . 5 12
12.0000 0.6783 0.3217 0.1095 6 11
13.0000 0.6166 0.3834 0.1156 7 10
14.0000* . . . 7 9
15.0000 0.5481 0.4519 0.1214 8 8
16.0000* . . . 8 7
17.0000 0.4698 0.5302 0.1268 9 6
18.0000* . . . 9 5
19.0000 0.3759 0.6241 0.1317 10 4
20.0000* . . . 10 3
21.0000 0.2506 0.7494 0.1348 11 2
24.0000 0.1253 0.8747 0.1113 12 1
32.0000 0 1.0000 0 13 0
* Censored Obrvation
Quantiles 75% 24.0000 Mean 17.1618
50% 17.0000 Standard Error 2.2053
25% 11.0000
这是用PL法对第2组生存资料的统计描述结果。各列的解释同上,从略。第2组患者的
中数生存期为17个月、平均生存期约为17个月。
Summary of the Number of Censored and Uncensored Values
GROUP Total Failed Censored �nsored
high-wbc 15 10 5 33.3333
low-wbc 20 13 7 35.0000
Total 35 23 12 34.2857
这是2组患者的总人数、死亡数 、删失数和删失百分比。
Survival Function Estimates
S SDF |
u 1.0 + *---HL
r | H*-HH-L---L
v | | L---L---L
i | H---H L-L-L
v | H-H L---L
a 0.5 + H-H L---L---L
l | H----H L---L
| H---H L-----L
D | | |
i | H------------H L---------------L
s 0.0 + H L
t ---+----+----+----+----+----+----+----+----+----+----+----+----+----+---
r 0.0 2.5 5.0 7.5 10.0 12.5 15.0 17.5 20.0 22.5 25.0 27.5 30.0 32.5
i
T
这是反映2组患者生存情况的生存曲线图,H表示有WBC倍增组、L表示无WBC倍增组,从图上可明显看出:无WBC倍增患者比有WBC倍增患者的生存期长。
四级能考几次 Censored Obrvations
Strata
相关表 L + L LL L L L L
H + H HH H H
-------+------+------+------+------+------+------+------+-------
0 5 10 15 20 25 30 35
T
这幅图反映了各组患者删失时间的情况。
-Log(Survival Function) Estimates
N -LOG SDF |
e 2 + H +L
g | + ++
a | + ++
t | + +L
i | +H ++
v 1 + ++ ++L
e | +H+ ++L+
| +H ++L+
L | H+++H ++L+L+L+
o | +H++L+++L+++L+
g 0 +*+++H*++
-+----+----+----+----+----+----+----+----+----+----+----+----+----+
S 0.0 2.5 5.0 7.5 10.0 12.5 15.0 17.5 20.0 22.5 25.0 27.5 30.0 32.5
D
T
这是按(t,-log(S(t))绘出的图,2条线都不呈直线趋势,说明生存时间不呈指数。
Log(-Log(Survival Function)) Estimates
L L(-L(S)) |
o 2 +
g |
| ++H +L
N 0 + +H++++H+ +L+L+L+
e | ++++H++H+ +LL++L+
g | +H+++ ++L+++L
a -2 + +++++H++++++L++++++L++
t | H+++*++++++
i |
v -4 +
e |
---+----+----+----+----+----+----+----+----+----+----+----+----+---
L 0.50 0.75 1.00 1.25 1.50 1.75 2.00 2.25 2.50 2.75 3.00 3.25 3.50
o
Log T
Legend for Strata Symbols
H:GROUP=high-wbc L:GROUP=low-wbc
这是按(logt, log(-log(S(t)))绘出的图,2条现别近似呈直线趋势,说明生存时间近似呈图尔。
Testing Homogeneity of Survival Curves over Strata
Rank Statistics
① GROUP Log-Rank Wilcoxon
high-wbc 5.073946 117
low-wbc -5.07395 -117
② Covariance Matrix for the Log-Rank Statistics
GROUP high-wbc low-wbc
high-wbc 3.36249 -3.36249
low-wbc -3.36249 3.36249
③ Covariance Matrix for the Wilcoxon Statistics
GROUP high-wbc low-wbc
high-wbc 2161.30 -2161.30
low-wbc -2161.30 2161.30
④ Test of Equality over Strata
Pr >
Test Chi-Square DF Chi-Square
Log-Rank 7.6565 1 0.0057
Wilcoxon 6.3337 1 0.0118
-2Log(LR) 2.8347 1 0.0922
这是关于各层生存曲线之间齐性检验的结果。①用2种检验方法算得公式(5.2.20)中的向量v,即v=(5.073946,-5.07395)'(Log-Rank法)、v=(117,-117)'(Wilcoxon法)。②、③分别用2种检验法算得此式中的协方差矩阵V,它们都是计算④中卡方值的中间结果,读者最需
要的是第④部分。这里给出了3种检验法的检验结果∶P值依次为0.0057(Log-Rank法)、0.0118(Wilcoxon法)、0.0922(似然比检验法,-2Log(LR))。
Life Table Survival Estimates
GROUP = high-wbc
⑴ ⑵ ⑶ ⑷ ⑸ ⑹
Conditional
Effective Conditional Probability
Interval Number Number Sample Probability Standard
[Lower, Upper) Failed Censored Size of Failure Error
0 5 4 1 14.5 0.2759 0.1174
5 10 3 4 8.0 0.3750 0.1712
10 15 2 0 3.0 0.6667 0.2722
15 20 1 0 1.0 1.0000 0
⑺ ⑻ ⑼ ⑽ ⑾
Survival Median Median
Interval Standard Residual Standard
[Lower, Upper) Survival Failure Error Lifetime Error
0 5 1.0000 0 0 9.1270 2.4177
5 10 0.7241 0.2759 0.1174 6.5000 2.1213
10 15 0.4526 0.5474 0.1440 3.7500 2.1651
15 20 0.1509 0.8491 0.1322 . .
Evaluated at the Midpoint of the Interval
⑿ ⒀ ⒁ ⒂
PDF Hazard
Interval Standard Standard
[Lower, Upper) PDF Error Hazard Error
0 5 0.0552 0.0235 0.064 0.031588
5 10 0.0543 0.0263 0.092308 0.051855
10 15 0.0603 0.0312 0.2 0.122474
15 20 0.0302 0.0264 0.4 0
这是用LT法对第1组资料进行统计描述的结果。编号⑴~⒂所代表的含义分别为:⑴按区间宽度=5将生存时间自动划分成若干区间、⑵死亡数、⑶删失数、⑷有效样本大小、⑸死亡的条件概率、⑹第⑸列数据的标准误差、⑺区间左端点处生存概率、⑻区间左端点处死亡概率、⑼第⑺列数据的标准误差、 ⑽中数剩余生存寿命(即在时刻ti活着的人有一半可望生存到的时间)、⑾第⑽列数据的标准误差、 ⑿区间中点概率密度函数的估计值、⒀第⑿列数据的标准误差、⒁区间左端点处危险概率的估计值、⒂第⒁列数据的标准误差。
此处是用LT法对第2组资料进行统计描述的结果(从略),解释方法同上。用3种检验方法给出的结果与用PL方法算得的结果相同,从略。
PL法与LT法的区别:在用计算机处理数据时,计算麻烦的苦恼已不存在,故PL法可适用于各种情况;当用手工计算且样本含量较大时,用LT法更方便一些。哩的计算结果基本上是一致的。PL法可看成是LT法的特殊情况,每个生存时间的区间宽度都为1。
[专业结论] 因前面的图示结果已表明此资料不服从指数,近似服从图尔,故宜选用Log-Rank法或Wilcoxon法检验的结果。均拒绝H侏儒角鲨0、接受H1,即2组生存曲线之间差别显著,
无WBC倍增患者的生存期显著地长于有WBC倍增的患者。
【例5.2.2】1965年某市肿瘤医院总结随访了15年曾在该医院住院手术的乳腺癌患者607例,整理后的资料如下,试分析该医院乳腺癌患者手术后的生存率。
术后年数: 0~ 1~ 2~ 3~ 4~ 5~ 6~ 7~ 8~ 9~ 10~
期内死亡人数: 59 69 43 30 13 7 14 4 3 0 0
小知不及大知 期内失访人数: 63 71 55 38 31 26 21 11 15 12 22
【分析与解答】此生存资料以分组的形式给出,SAS程序可按如下方式编写。
【SAS程序】━━【D5P2.PRG】
DATA abc; 63 59
KEEP f t c; 71 69
RETAIN t -0.5; 55 43
INPUT withdraw fail; 38 30
t=t+1; 31 13
c=0; 26 7
f=fail; 21 14
OUTPUT; 11 4
c=1; 15 3
f=withdraw; 12 0
OUTPUT; 22 0
CARDS; ;
(程序的第1部分) (程序的第2部分)
PROC LIFETEST PLOTS = (S, LS, LLS, H, P)
INTERVALS = ( 0 TO 10 ) METHOD = LT;
TIME t*c(1);
FREQ f;
RUN;
(程序的第3部分)
因篇幅所限,程序修改指导、输出结果及其解释等项内容从略。
第3节 生存资料非参数统计方法中的有关计算公式
1.乘积─极限法(Product-Limit,PL法或称为Kaplan-Meier,KM法)
让t1<t2<…<tk代表离散的失效(死亡或复发等)时间,设ni为第i个时刻开始之前生存的个体数目,即危险集的大小(i=1,2,…,k),再设di是在时刻ti失效的个体数目、si=ni-di。则在时刻ti的生存函数的PL估计值是ti时刻之前各时间点上生存率的乘积,即
(5.2.1)
式中的估计量属于右连续的,即在ti时刻发生的失效事件已包括在S^(ti)的估计中。与它对
应的标准误的估计值可用Greenwood的公式来计算:
(5.2.2)
生存时间的第1样本四分位数定义如下:
(5.2.3)
第2、3样本分位数可用类似的方式计算。q.50(即第2样本分位数)就是中位数,也就是中数生存期。
平均生存时间的估计值为∶
(5.2.4)
式中t0=0。若最后一个是删失数据,此式就低估了平均数。μ^的标准误差被定义为:
(5.2.5)
式中 。
2.寿命表法(Life-Table, LT法)
寿命表估计量通过计数落入时间区间[ti-1,ti]内的失效和删失的观察例数来计算,这
里i=1,2,…,k+1, t0=0、tk+1=∞。令ni为进入区间[ti-1,ti]内的个体数目、di和wi分别为发生在此区间内的事件(指死亡或失效,下同)数目和删失数目、bi=ti-ti-1、n'i=ni-wi/2;
n'i被称为此区间内有效的样本大小。再令tmi为此区间的中点、p^i=1-q^i。事件在此区间内发生的条件概率及其标准误差的估计值分别由式(5.2.6)和(5.2.7)定义:
(5.2.6) (5.2.7)
在时刻ti生存函数的估计值及其标准误差的估计值分别由式(5.2.8)和(5.2.9)定义:
(5.2.8)
(5.2.9)
在tmi处的密度函数及其标准误差的估计值分别由式(5.2.10)和(5.2.11)定义:
(5.2.10)
(5.2.11)
在tmi处的危险率函数及其标准误差的估计值分别由式(5.2.12)和(5.2.13)定义:
(5.2.12)
(5.2.13)
设在区间[tj-1,tj]内满足关系式S^(tj-1)≥S^(ti)/2>S^(tj),则在ti处的中数剩余寿命(即在ti活着的人有一半可望生存到的时间)及其标准误差的估计值由式(5.2.14)和(5.2.15)
定义:
微波炉多少瓦
(5.2.14)
(5.2.15)
3.各层间生存曲线的齐性检验
(1)对数秩检验和Wilcoxon检验
检验各层之间齐性所用的秩统计量为:
v'V-v~χ2 , df=R(V) (5.2.16)
式中R(V)(即协方差矩阵V的秩),从而可获得近似的概率水平。这里,v是一个c×1的向量(v1,v2,…,vc)',其具体表达式为:
(5.2.17)]
这里c是分层变量的层数;协方差矩阵V=(Vjl)由下式来估计:
(5.2.18)
这里下标i表示离散的失效时间,当j=l时,δjl=1;其他情况下,δ=0。nij和dij分别是第i个失效时间第j层中危险集的大小与事件数目, 、 、 。Vj可
以解释为:在生存曲线相同的假设之下,观察的与期望的失效数目之差的加权和。
当权wi=1时为对数秩检验;当权wi=ni时为Wilcoxon检验。V-代表V的广义逆矩阵。由此可知:数值较大的失效时间在对数秩检验统计量中所起的作用大;而数值较小的失效时间在威尔科克森检验统计量中所起的作用大(因通常寿命短的频数较大)。
(2)似然比检验
当资料各层之间服从指数的假设成立时,检验各层之间齐性(即检验指数的尺度参数相等)的似然比检验统计量由式(5.2.19)定义:
(5.2.19)
这里Nj是第j层内事件的总例数, , 是第j层中用于检验的总时间、mj是第j层中观察的总例数、 。把Z视为服从自由度为c-1的卡方分布,从而求得近似的概率水平。
4.协变量与生存时间联系密切程度的检验
用于检验协变量与生存时间联系密切程度的秩检验是用于齐性检验的秩检验的更一般推广,这种秩检验统计量具有如下的形式:
(5.2.20)
搞怪式中v由式(5.2.21)定义、V-为V的广义逆矩阵、V分别由式(5.2.22)与(5.2.24)定义:
(5.2.21)
当此式中 时,为对数秩得分检验;
当此式中 时,为Wilcoxon得分检验。
式(5.2.21)中下标及符号的含义如下:
α是观测对象的编号,α=1,2,…,n;n是观察的总数;i,j=1,2,…,k (k为不同的时间点数);t(j)代表有序的事件时间; Z(j)代表相应的协变量向量;tα代表有序时间(含删失与事件时间)。δα=1(如果观察到事件发生)、δα=0(如果观察到删失发生);得分Cα,δα取决于删失
的类型、且对全部观察值求和。
用于对数秩统计量的协方差矩阵的估计量是:
(5.2.22)
式中Vi是与t(j)时危险集所对应的校正的平和与交叉乘积和构成的矩阵,即
(5.2.23)
式中
用于Wilcoxon统计量的协方差矩阵的估计量为下式(其编号为∶(5.2.24))∶
此式中ai、a*i、Si、xi、si分别为:
、 、
、
第3章 COX模型回归分析
第1节 COX回归模型(半参数回归模型)
像通常的回归分析一样,人们也希望能建立起生存时间(因变量或反应变量)随危险因素(自变量或协变量)变化的回归方程,以便对危险因素的作用大小有一个全面的了解和掌握、并根据危险因素的不同取值对生存概率(或危险率)进行预测。由于生存时间的准确很难获得,前述目的很难直接实现。 1972年COX提出了比例危险模型(Proportional Hazard Model),简称为COX模型。由于此模型在表达形式上与参数模型相似,但在对模型中各参数进行估计时踊依赖于特定的假设,所以又有半参数模型之称。此模型的实用面
很宽,在生存分析中占有特殊的地位。其模型的具体形式为:
hi(t)=h0(t)exp(β1xi1+β2xi2+…+βmxim) (5.3.1)
式中hi(t)为第i名受试者生存到ti时刻的危险率函数,h0(t)是当所有危险因素(即xij=0)不存在时的基础危险率函数,X=(xi1,xi2,…,xim)'是可能与生存时间有关的m个危险因素所构成的向量。将式(5.3.1)变形如下:
ln[hi(t)/h0(t)]=β1xi1+β2xi2+…+βmxim (5.3.2)
此式表明:各危险因素与回归系数的线性组合就是第i名受试者的相对危险率函数的自然对数值。再设有i、j2个受试者,其危险因素向量分别为X1与X2,由式(5.3.1)不难得出他们的相对危险率的自然对数为:
ln[hi(t)/hj(t)]=β1(xi1-xj1)+β2(xi2-xj2)+…+βm(x中孝介花海im-xjm) (5.3.3)
即利用“具有某预后因素向量的受试者的死亡风险与不具有该预后因素向量的受试者的死亡风险在所有时间上都保持一个恒定比例”的假设,巧妙地获得了各时间点上2个受试者相对危险率函数的估计值。
然而,当资料不满足上述假设时,即有些危险因素作用的强度是随时间而变化的,2个受试者的危险率函数之比(相对危险)随时间而改变,就应改用时变协变量模型,也称为非比
例危险模型(Nonproportional hazard model)。当只有一个危险因素时,其模型的具体形式为:
hi(t)=h0(t)exp[βxi+γ(xiti)] (5.3.4)
式中ti为第i个受试者的生存时间。
上述各式中的回归系数需用最大似然法进行估计,一旦有了危险率函数的估计值,再利用生存时间函数之间的相互关系,可获得其他生存时间函数的估计值。
第2节 COX模型回归分析应用举例
[例5.3.1] 某医院肿瘤科提供的一份关于肺癌病人的失效时间资料,因变量(或反应变量)为病人治疗后的生存时间t(天),当t为删失数据时,用前面加一个负号来表示;考察的协变量(即危险或预后因素)如下。
①癌细胞的类型(Cell$),它有4个水平,即腺癌细胞(adeno)、鳞癌细胞(squamous)、小细胞肺癌(small)和大细胞肺癌(large);
②治疗类型(THERAPY$),它有2个水平,即标准的方法(standard)和试验的方法(test);
③疗前处理(PRIOR$),它有2个水平,即采取了疗前处理(yes)和未采取疗前处理(no);
④病人的年龄(age)(岁);
⑤从诊断到治疗的等待时间(diagtime);
ⅰ人的行动状态用Karnofsky率来度量, 其取值用KPS表示,10≤KPS≤30表明病人完全靠医院护理、40≤KPS≤60表明病人的行动部分地受到限制、70≤KPS≤90表明病人的行动可以自理。
前3个变量被当作分类变量,后3个变量被当作连续性变量。资料的形式为:
各组病人的治疗方法、癌细胞类型, 同一组中的样本含量
生存时间、KPS值、diagtime值、年龄、与疗前处理对应的指示变量PR值
( 注:PR=0等价于令PRIOR='YES',即表示采取了疗前处理; PR=10等价于令PRIOR='NO',即表示未采取疗前处理)。全部数据如下(文件名:lung.dat):
STANDARD SQUAMOUS 15
72 60 7 69 0 411 70 5 64 10 228 60 3 38 0 126 60 9 63 10
118 70 11 65 10 10 20 5 49 0 82 40 10 69 10 110 80 29 68 0
314 50 18 43 0 -100 70 6 70 0 42 60 4 81 0 8 40 58 63 10
144 30 4 63 0 -25 80 9 52 10 11 70 11 48 10
STANDARD SMALL 30
30 60 3 61 0 384 60 9 42 0 4 40 2 35 0 54 80 4 63 10
13 60 4 56 0 -123 40 3 55 0 -97 60 5 67 0 153 60 14 63 10
59 30 2 65 0 117 80 3 46 0 16 30 4 53 10 151 50 12 69 0
22 60 4 68 0 56 80 12 43 10 21 40 2 55 10 18 20 15 42 0
139 80 2 64 0 20 30 5 65 0 31 75 3 65 0 52 70 2 55 0
287 60 25 66 10 18 30 4 60 0 51 60 1 67 0 122 80 28 53 0
27 60 8 62 0 54 70 1 67 0 7 50 7 72 0 63 50 11 48 0
392 40 4 68 0 10 40 23 67 10
STANDARD ADENO 9
8 20 19 61 10 92 70 10 60 0 35 40 6 62 0 117 80 2 38 0
132 80 5 50 0 12 50 4 63 10 162 80 5 64 0 3 30 3 43 0
95 80 4 34 0
STANDARD LARGE 15
177 50 16 66 10 162 80 5 62 0 216 50 15 52 0 553 70 2 47 0
278 60 12 63 0 12 40 12 68 10 260 80 5 45 0 200 80 12 41 10
156 70 2 66 0 -182 90 2 62 0 143 90 8 60 0 105 80 11 66 0
103 80 5 38 0 250 70 8 53 10 100 60 13 37 10
TEST SQUAMOUS 20
999 90 12 54 10 112 80 6 60 0 -87 80 3 48 0 -231 50 8 52 10
242 50 1 70 0 991 70 7 50 10 111 70 3 62 0 1 20 21 65 10
587 60 3 58 0 389 90 2 62 0 33 30 6 64 0 25 20 36 63 0
357 70 13 58 0 467 90 2 64 0 201 80 28 52 10 1 50 7 35 0
30 70 11 63 0 44 60 13 70 10 283 90 2 51 0 15 50 13 40 10
TEST SMALL 18
25 30 2 69 0 -103 70 22 36 10 21 20 4 71 0 13 30 2 62 0
87 60 2 60 0 2 40 36 44 10 20 30 9 54 10 7 20 11 66 0
24 60 8 49 0 99 70 3 72 0 8 80 2 68 0 99 85 4 62 0
61 70 2 71 0 25 70 2 70 0 95 70 1 61 0 80 50 17 71 0
51 30 87 59 10 29 40 8 67 0
TEST ADENO 18
24 40 2 60 0 18 40 5 69 10 -83 99 3 57 0 31 80 3 39 0
51 60 5 62 0 90 60 22 50 10 52 60 3 43 0 73 60 3 70 0
8 50 5 66 0 36 70 8 61 0 48 10 4 81 0 7 40 4 58 0
140 70 3 63 0 186 90 3 60 0 84 80 4 62 10 19 50 10 42 0
45 40 3 69 0 80 40 4 63 0
TEST LARGE 12
52 60 4 45 0 164 70 15 68 10 19 30 4 39 10 53 60 12 66 0
15 30 5 63 0 43 60 11 49 10 340 80 10 64 10 133 75 1 65 0
111 60 5 64 0 231 70 18 67 10 378 80 4 65 0 49 30 3 37 0
[SAS程序]──[D5P3.PRG]
DATA valung; PROC PHREG;
RETAIN therapy cell; MODEL t*censor(1)=kps age diagtime;
LENGTH prior$ 3; RUN;
INFILE 'a:lung.dat'; PROC PHREG;
INPUT therapy$ cell$ n; MODEL t*censor(1)=kps age diagtime;
DO i=1 TO n; STRATA cell;
INPUT t kps diagtime age pr @@; RUN;
censor=(t<0); PROC PHREG;
t=ABS(t); MODEL t*censor(1)=kps age diagtime;
IF pr=10 THEN prior='yes'; STRATA cell therapy prior;
ELSE prior='no'; RUN;
OUTPUT;
END;
(程序的第1部分) (程序的第2部分)
[程序修改指导] 调用PHREG过程时,MODEL语句等号右边必须是连续性变量。这里所写的3个过程步的区别在于STRATA语句中所包含的枫变量的个数:第1个过程步不含STRATA语句,就是把所有资料看成来自1个总体;第2个过程步要求按CELL的4个水沏
分析资料;第3个过程步要求按CELL、THERAPY、PRIOR的16种水平组合枫分析资料。显然,在分层的条件下,COX所作的比例危险假设容易得到满足,但各层的样本含量不应太小。
另外,若连续变量很多时,可在MODEL语句最后增加选择项/SELECTION=方法名,进行变量筛选。方法名有如下几种:
BACKWARD或B(后退法)、FORWARD或F(前进法)、 STEPWISE或S(逐步回归法)、SCORE(最优回归子集法)。
[输出结果及其解释] The PHREG Procedure
Testing Global Null Hypothesis: BETA=0
Without With
Criterion Covariates Covariates Model Chi-Square
-2 LOG L 397.545 363.881 33.665 with 3 DF (p=0.0001)
Score . . 34.158 with 3 DF (p=0.0001)
Wald . . 31.510 with 3 DF (p=0.0001)
此结果表明:含3个自变量的COX回归模型是显著的(P<0.0001)。
Analysis of Maximum Likelihood Estimates
Parameter Standard Wald Pr > Risk
Variable DF Estimate Error Chi-Square Chi-Square Ratio
KPS 1 -0.036127 0.00655 30.44686 0.0001 0.965
AGE 1 -0.021571 0.01157 3.47602 0.0623 0.979
DIAGTIME 1 0.006426 0.01245 0.26624 0.6059 1.006
这是用最大似然法对模型中各参数估计并检验的结果,显然,只有变量KPS是显著的(P<0.0001)。这里给出的只是最后1个过程步输出的结果,因前2个过程步输出的结果与此结果相似,从略。详细的讨论请看下面的例子。
[例5.3.2] 若用LIFETEST过程预处理此资料会发现:不同细胞类型的生存曲线之间差别非常显著,生存曲线从左到右依次为:ADENO、SMALL、LARGE、SQUAMOUS。在[例5.3.1]中,请设法将变量CELL引入COX模型中来,重建模型。
[分析与解答] 因分类变量无法直接放入回归方程,这对模型的拟合是不利的。补救的办
法是:引入哑变量,使分类变量转变成数值袖量后再用PHREG过程(请看本例);若已知生存时间近似服从某特定的参数模型时,可直接用LIFEREG过程拟合参数模型,因参数模型中可包含2类变量(请看[例5.4.1])。于是,可对变量CELL作变换,使它变成3个哑变量。对数据文件作如下修改(改程序较困难):把ADENO改成1 0 0;把SMALL改成0 1 0;把LARGE改成0 0 1;把SQUAMOUS改成0 0 0。如第1行
STANDARD SQUAMOUS 15应改为:STANDARD 0 0 0 15,其他如法炮制。
程序[COXLUNG.PRG]的数据步中第1个INPUT语句和过程步需作如下修改
INPUT therapy$ adeno small large n;
程序[COXLUNG.PRG]中的过程步需作如下修改(因其他连续变量的作用不显著,故未将它们写入下面的模型语句之中):
PROC PHREG;
MODEL t*censor(1)=kps adeno small large;
STRATA therapy prior;
RUN;
[输出结果及其解释]
The PHREG Procedure
Testing Global Null Hypothesis: BETA=0
Without With
Criterion Covariates Covariates Model Chi-Square
-2 LOG L 694.500 641.270 53.230 with 4 DF (p=0.0001)
Score . . 55.507 with 4 DF (p=0.0001)
Wald . . 52.163 with 4 DF (p=0.0001)
Analysis of Maximum Likelihood Estimates
Parameter Standard Wald Pr > Risk
Variable DF Estimate Error Chi-Square Chi-Square Ratio
KPS 1 -0.031260 0.00548 32.50520 0.0001 0.969
ADENO 1 1.106441 0.30933 12.79392 0.0003 3.024
SMALL 1 0.806382 0.27902 8.35214 0.0039 2.240
LARGE 1 0.364903 0.29681 1.51148 0.2189 1.440
关于模型中全部参数的检验结果是非常显著的(P<0.0001)。变量KPS的作用是非常显著的(χ2=32.505,P<0.0001),其相对危险度为0.969,即KPS取得1个正的增量时,相对危险度会有所降低;与SQUAMOUS类型的癌细胞相比:ADENO和SMALL都非常显著,相对危险度分别为3.024和2.240倍;LARGE不显著,其相对危险度是SQUAMOUS的1.44倍。
[专业结论] 说明病人的行动状态KPS取值大,生存时间长;生存时间由短到长的癌细胞类型依次为(根据它们的回归参数和相对危险度的估计值判断):ADENO、SMALL、LARGE、SQUAMOUS。
第4章 参数模型回归分析
第1节 参数回归模型
通过图解法或统计检验法,若能得知欲分析的生存资料服从某特定的参数模型,如指数模型或图尔模型时,直接拟合这些模型,便可获得比盲目用其他法更准确的结果。现将几种
最常用的参数模型所对应的生存时间函数的具体形式汇总如下:
参数模型名称 生存概率 密度函数 危险函数
S(t) f(t) h(t)
━━━━━━ ━━━━ ━━━━━━━━ ━━━━━━━
指数分布 ρ
威布尔分布