SAS第十课：随访资料的生存分析－－非参数法与半参数Cox比例风险模型
2012-04-17 生物谷 生物谷
生存分析方法大体上可分为三类：非参数法、参数法和半参数方法，与之相对应，SAS提供了三个程序步用于生存分析，它们是： LIFETEST过程 提供非参数分析方法，用乘积极限法(Product limit method)和寿命表法(Life table method)估计生存率和中位生存时间等；用对数秩检验(Log-rank test)、Wilcoxon检验和似然比检验等做分组
生存分析方法大体上可分为三类：非参数法、参数法和半参数方法，与之相对应，SAS提供了三个程序步用于生存分析，它们是：
- LIFETEST 过程 提供非参数分析方法，用乘积极限法(Product limit method)和寿命表法(Life table method)估计生存率和中位生存时间等；用对数秩检验(Log-rank test)、Wilcoxon检验和似然比检验等做分组比较。该过程主要用于估计生存率及进行单因素分析。
- LIFEREG 过程 提供指数模型、Weibull模型、Gompertz模型等参数分析方法。
- PHREG 过程 提供半参数Cox比例风险模型分析。
本章只介绍常用的LIFETEST过程和PHREG过程。
§10.1 LIFETEST过程及其应用
10.1.1 语法格式
LIFETEST
过程的语法格式如下：|
PROC LIFETEST [选项]; |
|
TIME <生存时间变量*截尾指示变量(数值)> ; |
|
[TEST <分组变量名列>; |
|
STRATA <分组变量名列>; |
|
FREQ <变量名列>; |
|
BY <变量名列>;] |
10.1.2 语法说明
PROC
和TIME语句为必需的，其他语句都可以省略。TIME语句，为必需语句,定义生存时间和截尾指示变量【过程选项】
- PL
- PLOTS=
【
TIME语句】TIME
语句用于定义生存时间和截尾指示变量。对截尾指示变量可以指定发生失效事件的数值，默认失效事件用0来表示，截尾事件用1来表示。【
STRATA语句和TEST语句】STRATA
语句定义生存率比较的分组变量，TEST语句定义生存率比较的分组变量或协变量。STRATA语句在这里的作用和BY语句类似，都是要求按分组变量名列进行分析，在计算生存率时各组分开计算；TEST
语句定义需检验的变量，即生存时间与该变量是否有关，如果它后面定义的变量为数值变量，则把该变量当作协变量检验与生存时间的关系。如果它定义的为分组变量，则分组比较生存时间有无差别。10.1.3 应用实例
例
10.1 观察两组卵巢腺癌患者的病程天数如下。请用对数秩检验比较两组的生存期差异有无统计学意义，并作生存率曲线。(医统P.321，7.2题)A
组(低恶性高分化癌)：28 29 175 195 309 377＋ 393＋ 421＋ 447＋ 452 709＋ 744＋ 770＋ 1106＋ 1206B
组(高恶性低分化癌)：34 88 137 199 280 291 299＋ 300＋ 309 351 358 369 370 375 382 392 429＋ 451 1119＋解：程序如下：
|
data a.yt7_2; |
|
input t @@; |
|
if t<0 then censor=1; |
|
else censor=0; |
|
if _N_<16 then group='A'; |
|
else group='B'; |
|
t=abs(t); |
|
cards; |
|
28 29 175 195 309 -377 ... |
|
... 382 392 -429 451 -1119 |
|
; |
|
proc lifetest method=pl plots=(s,ls); |
|
time t * censor(1); |
|
strata group; |
|
proc lifetest method=lt plots=(s,h); |
|
time t * censor(1); |
|
strata group; |
|
run; |
本例为不分组资料，生存时间有截尾值，应用乘积极限法和寿命表法分别作生存分析，结局变量只有一个，即生存时间，负值表示截尾，用
if语句完成转换，产生新变量CENSOR，其取值1表示截尾，0表示不截尾。用两个PROC语句，分别用不同的方法作生存分析，按GROUP变量分组统计。LIFETEST
过程的默认输出结果很多，主要有：① 生存率及其标准误。
② 生存时间四分位数。
③ 生存时间的均值及其标准误。
④
Wilcoxconc2统计量。当生存时间符合对数正态分布时该检验统计效率较高。⑤
Log-rankc2统计量。当生存时间符合Weibull分布或属于比例风险模型时该检验统计效率较高。本例程序运行的主要结果如下：
第一部分，PL法对生存资料分组进行统计描述
The LIFETEST Procedure
Product-Limit Survival Estimates
GROUP = A
Survival
Standard Number Number
T Survival Failure Error Failed Left
生存时间 累积生存率 死亡概率 累积生存率标准误 已观察到的不同 观测尚未结束的例数
失效时间的例数
0.00 1.0000 0 0 0 8
28.00 0.9333 0.0667 0.0644 1 7
29.00 0.8667 0.1333 0.0878 2 6
175.00 0.8000 0.2000 0.1033 3 5
195.00 0.7333 0.2667 0.1142 4 4
309.00 0.6667 0.3333 0.1217 5 3
377.00* . . . 5 2
393.00* . . . 5 1
1206.00 0 1.0000 0 6 0
* Censored Observation 标有"*"的为截尾值 Summary Statistics for Time Variable 时间变量名
生存时间四分位数、点估计及95%可信区间
Point 95% Confidence Interval
Quantile Estimate [Lower, Upper)
75% 1206.00 195.00 1206.00
50% 252.00 34.00 1206.00
25% 34.00 28.00 309.00
Mean 458.87 Standard Error 164.95
生存时间均数 均数的标准误
Summary of the Number of Censored and Uncensored Values
各组患者的总人数、死亡数、截尾数和截尾百分比
GROUP Total Failed Censored %Censored
分组变量 总人数 死亡数 截尾数 截尾百分比
A 10 8 2 20.0000
B 14 14 0 0.0000
Total 24 22 2 8.3333 Testing Homogeneity of Survival Curves over Strata
各组所在总体生存时间是否相等的检验
Time Variable 时间变量
Test of Equality over Strata
Pr >
Test Chi-Square DF Chi-Square
统计量 值 自由度 P值
Log-Rank 5.1240 1 0.0236
Wilcoxon 1.9548 1 0.1621
-2Log(LR) 4.0098 1 0.0452
这是各层生存曲线之间齐性检验的结果。这里给出了3种检验方法的结果：P值依次为0.0236(log-rank test)、0.1621(Wilcoxon test)、0.0452(likelihood ratio test)。此外，还给出log-rank test和Wilcoxon test检验中卡方值的中间结果(本处未列出)。按α=0.05水准，拒绝H0，接受H1 ，两组生存时间有差别。
第二部分是寿命表法结果，把生存时间划分为区间，计算区间左端点处的生存概率。
Life Table Survival Estimates
GROUP = A
Conditional
Effective Conditional Probability
Interval Number Number Sample Probability Standard
[Lower, Upper) Failed Censored Size of Failure Error
(1) (2) (3) (4) (5) (6)
按区间宽度为200将 死亡数 截尾数 校正期初人数 死亡的条件概率 第(5)列数据
生存时间划分成若干区间 的标准误 0 200 4 0 15.0 0.2667 0.1142
200 400 1 2 10.0 0.1000 0.0949
400 600 1 2 7.0 0.1429 0.1323
600 800 0 3 3.5 0 0
800 1000 0 0 2.0 0 0
1000 1200 0 1 1.5 0 0
1200 1400 1 0 1.0 1.0000 0 Survival Median Median
Interval Standard Residual Standard
[Lower, Upper) Survival Failure Error Lifetime Error
(1) (7) (8) (9) (10) (11)
按区间宽度为200将 区间左端点处 区间左端点处 第(7)列数据 中位剩余生存时间 第(10)列数据
生存时间划分成若干区间 生存率 死亡概率 的标准误 的标准误
0 200 1.0000 0 0
200 400 0.7333 0.2667 0.1142
400 600 0.6600 0.3400 0.1241
600 800 0.5657 0.4343 0.1376
800 1000 0.5657 0.4343 0.1376
1000 1200 0.5657 0.4343 0.1376
1200 1400 0.5657 0.4343 0.1376Evaluated at the Midpoint of the Interval
PDF Hazard
Interval Standard Standard
[Lower, Upper) PDF Error Hazard Error
(1) (12) (13) (14) (15)
区间中点概率密度 第(12)列数据的标准误 区间左端点处风 第(14)列数据
函数的估计值 险率估计值 的标准误
0 200 0.00133 0.000571 0.001538 0.00076
200 400 0.000367 0.000353 0.000526 0.000526
400 600 0.000471 0.000445 0.000769 0.000767
600 800 0 . 0 .
800 1000 0 . 0 .
1000 1200 0 . 0 .
1200 1400 0.00283 0.000688 0.01 0
给出的其它结果同
log-rank test，此处从略。两法计算结果基本一致，PL法可看成是LT法的特殊情况，即每个生存时间的区间宽度为1。§10.2 PHREG过程及其应用
PHREG
（proportional harzard regression,比例风险回归）过程基于Cox比例危险模型对生存数据进行回归分析，结局变量（应变量）为生存时间，可以处理生存时间有截尾的数据。模型中的自变量可以是连续性、分类变量、时间依存的自变量。对成比例风险是否成立作出检验，利用最大似然法迭代求出模型的参数估计，对模型的参数作似然比、比分检验和Wald检验三种检验。10.2.1 语法格式
PHREG
过程的语法格式如下：|
PROC PRREG [选项]; |
|
MODEL <生存时间变量*截尾指示变量(数值)>=<自变量名> /[选项];[1] |
|
[STRATA <分组变量名列>; |
|
FREQ <变量名列>; |
|
BY <变量名列>;] |
10.2.2 语法说明
1
．MODEL语句为必需的,定义生存时间和截尾指示变量和说明变量【过程选项】
【模型选项】
- TIES= 方法 指定估计生存率所用的方法：
- BRESLOW 使用Breslow的近似似然估计，为默认的选项
- DISCRETE 用离散Logistic模型替代比例风险模型，多用于m:n的Logistic回归
- EFRON 使用Efron的近似似然。
- EXACT 计算在比例危险假定下所有失效事件发生在具有相同值的删失时间或较大值时间之前的精确条件概率。
- ENTRYTIME= 变量名，规定一个替代左截断时间的变量名。
- SELECTION=method, 方法可以选择以下几种，
- FORWARD( 或F) 按照规定的P值SLE从无到有依次选一个变量进入模型
- BACKWARD （或B） 按照规定的P值SLS从含有全部变量的模型开始，依次剔除一个变量
- STEPWISE （或S） 按照SLE的标准依次选入变量，同时对模型中现有的变量按SLS的标准剔除不显著的变量
- SCORE 采用最优子集选择法
【STRATA语句】
比例风险的假定可能不会对所有的层都成立，此时需要作分层分析。STRATA语句要求按照分层变量名列的水平数拟合一个多层的Cox模型。与BY语句不同，后者是要求按分组变量名列分别估计模型及参数。
PHREG
过程中还可以加入编程语句用以创建模型中的新的自变量，但不能用以修改应变量，截尾变量，组变量或分层变量的值。当省略所有的选项，并且只有一个分类自变量（分组变量）时，模型的检验相当于生存曲线的比较。10.2.3 应用实例
例10.2 随访25例分别以A、B治疗方法治疗的某癌症病人，资料如下，+号表示为截尾值。1：有肾功能损害，0：无肾功能损害，请试作COX回归。
|
A疗法 |
B疗法 | ||||
|
编号 |
肾功损害 |
生存日数 |
编号 |
肾功损害 |
生存日数 |
|
1 |
1 |
8 |
13 |
1 |
13 |
|
12 |
0 |
52 |
16 |
1 |
18 |
|
5 |
1 |
58 |
25 |
1 |
23 |
|
8 |
1 |
63 |
11 |
0 |
70 |
|
21 |
1 |
63 |
10 |
0 |
76 |
|
7 |
0 |
220 |
2 |
0 |
180 |
|
24 |
0 |
365 |
9 |
0 |
195 |
|
4 |
0 |
452 |
20 |
0 |
210 |
|
18 |
0 |
496 |
3 |
0 |
232 |
|
22 |
0 |
528+ |
17 |
0 |
300 |
|
19 |
0 |
560+ |
23 |
0 |
396 |
|
15 |
0 |
676+ |
14 |
0 |
490+ |
|
6 |
0 |
540+ | |||
利用PHREG过程拟合COX比例风险模型，建立一个数据集，生存日数day为结果变量，还需一个截尾指示变量censor，1为截尾，0为无截尾。组变量GROUP指示治疗组，协变量renal表示肾功能损害。
程序如下：
|
data a.bk5_2; |
|
input group renal day; |
|
censor=(day<0); |
|
days=abs(day); |
|
cards; |
|
1 1 8 |
|
... |
|
2 0 -540 |
|
； |
|
proc phreg data=a.bk5_2; |
|
model days*censor(1)=group renal ; |
|
run; |
定义截尾指示变量censor，用逻辑运算函数实现转换，当day<0为真时，取值为1，否取值0，对day取绝对值为了使参加运算的数都是正值。调用PHREG过程，以
days为应变量，自变量为处理组和肾功能损害拟合COX模型。程序运行的主要结果如下：The PHREG Procedure
Data Set: A.BK5_2 数据集名称
Dependent Variable: DAYS 应变量名
Censoring Variable: CENSOR 截尾指示变量
Censoring Value(s): 1 截尾值
Ties Handling: BRESLOW BRESLOW法处理相等的数据 Summary of the Number of
Event and Censored Values Percent
Total Event Censored Censored 25 20 5 20.00
总例数 死亡数 截尾数 截尾的百分数 Testing Global Null Hypothesis: BETA=0
模型检验，无效假设为β=0
， Without With ， ，
Criterion Covariates Covariates Model Chi-Square
-2 LOG L 106.176 83.260 22.916 with 2 DF (p=0.0001)似然比检验
Score . . 29.715 with 2 DF (p=0.0001)比分检验
Wald . . 13.863 with 2 DF (p=0.0010) Wald检验 Analysis of Maximum Likelihood Estimates
参数的最大似然估计
， ， Parameter Standard Wald Pr> Risk Variable DF Estimate Error Chi-Square Chi-Square Ratio 变量名 自由度 参数估计 标准误 参数的Waldχ2检验 P值 相对危险度 GROUP 1 0.989726 0.52355 3.57363 0.0587 2.690 RENAL 1 4.112210 1.13854 13.04529 0.0003 61.082
本例模型总的检验三种方法的P值都小于0.05，模型有统计学意义。
对自变量的检验结果用Waldχ2检验，P值分别为0.0587，0.0003。根据参数估计值可写出COX回归方程：
h(t,x)=h0(t)*e0.989726group+4.112210renal相对危险度分别为2.690，61.082，说明B组死亡的危险为A组的2.690倍，而伴肾功能损害的死亡的危险为无肾功能损害61.082倍。
（原著：田晓燕 李晓松）
本网站所有内容来源注明为“梅斯医学”或“MedSci原创”的文字、图片和音视频资料，版权均属于梅斯医学所有。非经授权，任何媒体、网站或个人不得转载，授权转载时须注明来源为“梅斯医学”。其它来源的文章系转载文章，或“梅斯号”自媒体发布的文章，仅系出于传递更多信息之目的，本站仅负责审核内容合规，其内容不代表本站立场，本站不负责内容的准确性和版权。如果存在侵权、或不希望被转载的媒体或个人可与我们联系，我们将立即进行删除处理。
在此留言
#非参数#
76
#SAS#
67
#风险模型#
66