分段指数模型的拟合及SAS实现

562

中国卫生统计2009年12月第26卷第6期

分段指数模型的拟合及SAS实现

山西医科大学卫生统计教研室(030001)

杨珊珊 彭江滨 余红梅 陈 芊

提 要 目的 以实例阐述生存分析分段指数模型的拟合及SAS实现。方法 将生存时间划分为几个区间,利用SAS中的PROCGENMOD过程或PROCLIFEREG过程,对生存资料进行分段指数模型拟合,并通过参数估计结果,计算各时间区间的风险率及生存率。结果 该资料拟合分段指数分布的结果比较满意,三个时间区间的死亡风险均不相同,且第一个和第二、三个区间的差别有统计学意义(P<0 05)。结论 分段指数模型不仅模型形式简单,参数易于估计,若拟合得当可提高统计效率。

关键词 生存分析 分段指数模型 生存率 风险率

指数分布是历史上第一个生存时间分布,是最简单与最重要的分布,在生存分析中占据重要地位。其特点为在任何时间上的风险函数为一常数,但实际数据通常不满足此要求。如果生存资料满足两个假设:每个时间区间的失访是均匀分布的,而且与治疗失败无关;每个时间区间内发病或死亡的概率较小,并且每个区间的发病或死亡时间服从独立的指数分布,则可以拟合分段指数模型。本文将结合实例介绍分段指数模型在生存分析中的应用及其SAS实现。

原理与方法

将生存数据的时间尺度划分为几个区间,假设各区间内的风险是常数,但区间之间的风险可不同。用符号表示如下,定义j个区间,其分割点为a1,a2, ,aj,a0=0,aj= ,区间j为[aj-1,aj)。分段指数模型假设每个时间区间有独立的指数分布,其个体i的风险函数如下:

hi(t)= je

xi

1

系: j=- j

2

。

实例分析与SAS实现

资料来源于17年里对149位糖尿病患者随访的

3

临床病例,变量及其赋值如下:

表1 糖尿病资料中变量及其编码

变量名称结局生存时间随访开始时年龄体重指数

诊断出糖尿病时的年龄吸烟状况收缩压舒张压心电图读数病人是否有冠心病

代码deadstimeAge1BMIAge0smkSBPDBPECGCHD

岁

不吸烟=0,曾吸烟=1,吸烟=2mmHgmmHg

正常=0,可疑=1,异常=2无=0,有=1编码

死亡=1,截尾=0年岁

*:考虑到收缩压和舒张压两个变量有一定的相关性,数据分析时取平均血压(MBP),MBP=SBP (1/3)+DBP (2/3)。

先对该资料拟合指数回归、Weibull回归、标准

(1)(2)

Gamma回归和Gamma回归模型,结果表明七个变量中与糖尿病生存有关的变量为随访开始时年龄、平均血压和心电图读数,且拟合优度检验拒绝指数回归模型,所以考虑拟合分段指数模型,只纳入前述三个协变量及指示变量j。

将糖尿病患者的生存时间划分为三个区间,生成一个新数据集year,然后用广义线性模型的GENMOD过程估计分段指数模型。GENMOD过程的优点是可以在MODEL语句中自由指定交互作用。在GEN-MOD过程中,需要在DATA步中先定义一个对数时间变量,本例在OUTPUT语句前加入ltime=log(time)。因为分段指数模型的最大似然估计值与Poisson回归的最大似然估计值相同,所以MODEL语句的分布选择项中选用poisson分布(dist=poisson)。link函数选用对数函数,offset规定了产生死亡危险集合的是对数时间变量。GENMOD过程报告的是应变量为对数风险函数的结果。具体SAS过程如表2所

4

其中aj-1 t<aj,或等价: loghi(t)= xij+

其中 j=log j。因此,允许对数风险方程的截距在区间之间变化,且这种变化没有固定的形式。

根据模型的参数估计结果可进一步计算生存率,各时间区间生存率的计算方法如下

- jbj

1

:

(3)

第j时间区间内的生存率=第(j-1)时间区间的生存率 e

其中bj为第j个时间区间的长度。

SAS软件中的PROCGENMOD过程或PROCLIFEREG过程可实现分段指数模型的估计。PROCGENMOD过程报告的是应变量为对数风险函数的结果,PROCLIFEREG过程报告的是应变量为对数生存时间的结果,两结果完全等价,且系数有以下关

通讯作者:余红梅

1,2

ChineseJournalofHealthStatistics,Dec2009,Vo.l26,No.6

563

示,拟合结果见表3。

表2 糖尿病数据的分段指数模型及SAS程序

(GENMOD过程)

行号010203040506070809101112131415161718

程序

DATAyear;SETdiabetes;

year=ceil(stime/6);DOj=1TOyear;time=6;event=0;IFj=yearanddead=1THENDO;event=1;time=stime-6*(year-1);END;ltime=log(time);OUTPUT;END;RUN;

PROCGENMODdata=year;CLASSjECG;

MODELevent=AgelECGMBPj/dist=poissonlink=logoffset=ltime;RUN;

S(time)=exp[-time exp(-11 3402+0 0865x1-1 9671x2+0 0452x3-2 0125J1)]其中x1为变量AGE1,x2为变量ECG0(心电图读数正常),x3为变量MBP,J1为生存时间三个区间的第一个区间。

由表3结果可知指示变量j(年)对应的P=0 0036<0 05,有统计学意义,说明风险函数是随时间而变化的,但在具体的某一时间区间内,其风险函数仍为一常数。由于j1和j2的系数都是以j3为参照组,因此可得出第一个6年的死亡风险明显低于第三个6年,且差别有统计学意义;第二个6年的死亡风险低于第三个6年,但两个区间风险的差别没有统计学意义。j1与j2比较的Wald检验由协方差矩阵构建,公式为

=

2

2

:(4)

( j- j)12

2

Var( j1)+Var( COV( j2)-2j1j2)

其中 j1, j2分别为j1,j2的系数估计值,Var( j1),Var( j2)分别为j1,j2自身相交处的方差值,COV( j1 j2)为j1,j2相交处的协方差值。带入相应数值得 =7 9680, =1,得P<0 05,故第一个6年与第二个6

年的死亡风险差别有统计学意义,第一个6年的死亡风险低于第二个6年。

2

由表3可得到相应的风险函数为:

=exp(-11 3402+0 0865x1-1 9671x2+0 0452x3-2 0125J1)

相应的生存函数为:

表3 糖尿病资料分段指数模型拟合结果

变量InterceptAgelECGECGECGECGjj1j2j3ScaleWeibullShape

012

自由度1121101211000

系数-11.34020.0865 -1.9671-1.1948-0.00000.0452 -2.0125-0.5780-0.00001.00001.0000

标准误

下限

2.16690.0196 0.56080.64070.00000.0168 0.6454

0.54180.00000.00000.0000

-15.58720.0482 -3.0662-2.45050.00000.0123 -3.2775

-1.63990.00001.00001.0000

95%CI

上限-7.09320.1249 -0.86810.06100.00000.0780 -0.7675

0.48400.00001.00001.0000

Wald 227.3919.5512.3912.313.48 7.2411.269.72

1.14

P值<0.0001<0.00010.00200.00050.0622 0.00710.00360.0018

0.2861

MBP

表4 分段指数模型估计的生存率

心电图读

时间区间

数分组正常正常正常可疑可疑可疑异常异常异常

[0,6)[6,12)[12,18)[0,6)[6,12)[12,18)[0,6)[6,12)[12,18)

风险率公式e^ +^ 1+^j1e^ +^ 1+^j2

^+^ 1

e

风险率估计值生存率估计值

0.00150.00640.01130.00330.01380.02450.01080.04540.0810

0.99090.96260.93430.98050.92080.86320.93710.76140.6151

Age1为50岁(随访开始时年龄中位数) …… 此处隐藏：5353字，全部文档内容请下载后查看。喜欢就下载吧 ……