11生存分析与SAS程序

在医学研究中，考察研究因素对研究对象的效应，经典的研究设计方法只观察其所出现的结局或结果。但是在一些诸如癌症等险恶疾病和慢性病的随访研究中，只观察其结局就不够全面，还需要考察研究对象出现某种结局所经历的时间，这段时间称为生存时间（survivaltime）。生存分析是用来分析生存时间资料的统计方法，是近几十年来产生且发展甚为迅速的一门应用统计的分支。到目前为止，生存分析已形成一套完整的体系，包括描述生存规律和进行组间比较的参数和非参数方法以及分析影响生存期因素的回归模型等。本章主要介绍描述生存规律的寿命表和乘积限估计法，用于组间比较的logrank检验和作预后因素分析的Cox回归模型以及相应的SAS过程。11.1生存分析的基本概念11.1.1生存时间在医学研究中，对于肿瘤、心血管等慢性疾病，要考察其治疗方法优劣，疾病预后的好坏以及影响疾病预后的因素，通常采用随访研究的方法。对某一疾病作随访研究时，一般是从某一时间开始，观察到某一规定时间截止，而研究对象是始点以后陆续进入观察。随访中要规定一个事件作为随访结局，例如：病人死于研究疾病。如果病人的随访结果是规定的结局，则称为失效或死亡（Failure），那么病人从进入观察随访到规定的结局出现，其间所经历的这段时间称为生存时间或失效时间（Survivaltime,Failuretime,Waitingtime）。通常用Ti表示（i=1,2,n,为观察个体数）。在随访中，由于客观条件限制，不能也不可能将全部观察对象都观察到规定的结局（Failure）。在观察截止时，一组观察对象除了出现规定的结局外，还有三种结果：（1）死于其他疾病；（2）由于迁移等原因失去联系（失访）；（3）随访截止时尚未出现规定的结局。这三种结果虽原因不一，但提供的信息是一致的，这类病人不能获得确切的生存时间，而只知道其生存期比随访观察到的时间长的信息。我们将此类病例数据称为截尾（Censored）数据。其观察到的时间用G表示。在获得生存时间时要注意以下几点：随访结局可以不是死亡。随访的结局根据研究问题的性质和目的的不同可以是死亡、复发、恢复等。但要明确时间界线。例如：肿瘤病人已作了手术切除治疗，我们随访观察他有否复发，复发就是我们规定的结局。1. 确定进入随访观察的起点时间。随访研究中，无论是同时进入观察还是陆续进入观察，对所有对象都要规定某一事件作为进入观察起点事件。例如：病人确诊，治疗开始或手术时间，动物试验中染毒等，作为起点事件要有明确的时间界线。3时间尺度。医学研究一般用日历时间作为时间尺度。但在某些医学问题研究中，可能不以日历时间作为时间尺度更能反映所研究的问题。例如：对儿童龋齿出现的随访研究中常以儿童出到第几颗牙才出现龋齿作为时间尺度。11.1.2生存数据的特征生存数据在结构上有它自己的特点，主要表现为：生存时间数据的分布与常见的统计数据分布有明显不同，常呈指数分布、Weibull分布、对数正态分布、对数Logistic分布、Gamma分布或更为复杂的分布，生物医学中的生存时间分布有时呈现不规则状态，因而难以用传统的统计方法对这类数据进行处理。1. 生存数据中往往包含有截尾数据。在随访研究中，若能观察到研究所规定的结局就能获得确切的生存时间，这类数据称为完全数据。但在实际工作中，由于时间和客观条件限制，部分患者难以观察到终点事件，以致不能获得确切的生存时间，而只能获得进入观察至失访时这段时间。在分析总结时若剔除这部分对象往往导致样本有偏性，并且这部分数据仍然提供了其实际生存时间大于观察到的截尾时间的信息，应充分利用这部分资料的信息。因此，一般随访研究的样本数据中不可避免地包含有截尾数据。11.1.3生存时间函数描述生存时间分布规律的函数主要有生存函数、死亡概率函数、概率密度函数和危险率函数。为了后文叙述方便，这里主要介绍生存函数和危险率函数。生存函数在描述生存规律的数量指标中，以往常用的指标是某个特定时间的生存率(例如：3年生存率、5年生存率)。这一指标的主要缺陷为不能反映整个生存规律，一个理想的指标应该是任意时间的生存率，即生存率是任意时刻t的函数。其意义是研究个体生存时间长于t的概率。若令T为生存期，s(t)为任意时刻t的生存率，得s(t)=p(Tt)0wtg(11.1)则称s(t)为生存率函数，简称生存函数。(SurvivorFunction)。将s(t)对t在直角坐标系作图画出生存率曲线，可从图上粗略估计出中位生存期或半数生存期。即生存曲线上取生存率为50%时所对应的时间。危险率函数(HazardsFunction)如果我们考虑已活到t时刻的患者，在时间t附近的瞬间死亡危险性，根据数学上极限性质，可表示为：lim(t兰Tct+心t|T3t)h(t)=.(11.2)心t0At则称h(t)为危险率函数。即相当于条件瞬间死亡率或年龄别死亡率。这是生存分析中一个特有的函数，为处理生存分析问题带来了极大方便。s(t)与h(t)的关系h(u)du亠上，亠/s(t)=e0=e其中H(t)=oh(u)du，称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。首先计算第一个时间区间to-的n。、do和Wo，然后计算下一个区间的门仆片=no-do-w。，以此类推，得：ni+i=ni-di-Wi(11.3)表11.1296例晚期肝癌患者生存率计算(时间单位：月)随访月数期初观期内死期内失访校正人死亡概生存概生存率生存率的标准ti(1)察人数亡人数人数Wi数ni率qi率Pis(ti)(8)误SEs(ti)ni(2)di(3)(4)(5)(6)(7)(9)0-29694102910.32300.67700.67700.02741-1927415184.50.40110.59890.40550.02942-1032210980.22450.77550.31440.02853-71226680.32350.67650.21270.02634-435540.50.12350.87650.18640.02555-3366300.20000.80000.14920.02466-214120.50.19510.80490.12010.02377-162115.50.12900.87100.10460.02308-1332120.25000.75000.07850.02179-82080.25000.75000.05890.020210-62250.40000.60000.03530.017711及以上2221.00000.00000.0000(3)校正观察人数ni的计算由于在ti-区间内的失访人数Wi在该区间内并未观察至区间的终点，因而该区间内的有效人数并非为ni。设Wi个失访个体在该区间内为均匀分布，即假定Wi个体平均观察半个区间，那末校正人数为：mnWi/2(11.4)(4) 死亡概率q和生存概率P计算qi是指在ti-区间始端仍活着的病人在该区间内死亡的可能性，其计算公式为：qi=di/n/(11.5)本例qi计算见表11.1第(6)栏。pi是指在ti-区间始端仍活着的病人在该区间内生存的概率，其计算公式为：Pi=1-qi(11.6)本例计算结果见表11.1第(7)栏。(5) 生存率s(ti)计算。根据概率的乘法法则，s(ti)的计算公式为：s(ti)=PoP1P2Pi-1(11.7)本例计算结果见表11.1第(8)栏(6) 生存率标准误SEs(ti)的计算。生存率是依据实际资料(样本)计算得到的，故有抽样波动，须对其标准误SEs(ti)进行估计，其近似计算公式为：SEs(ti)=s(tjqo_“_q1_.“qy.(11.8)VonoP1Pgn求得生存率和标准误，就可以对这批随访对象在整个随访期间的生存规律进行描述；并可按正态近似原理，对总体生存率的可信区间作出估计。乘积限估计法当数据个体较少时，为充分利用每个数据的信息，必须采用更为精确的估计方法。这些估计方法中应用最多、效率较高的是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)ti（2）山（3）di（4）qi（5）Pi（6）s(t)(7)SEs(ti)(8)12160011241510.06670.93330.933336140010.93330.064446130010.933357.5120010.933368.51110.09090.90910.84850.099979100010.8485810910.11110.88890.757291280010.75420.12571013710.14290.85710.64640.14681118610.16670.83330.53870.1569121950010.53871324410.25000.75000.40400.16571426310.33330.66670.26930.15581531210.50000.50000.13470.123116431001（2）计算观察人数口和死亡人数dj。先分别列出各死亡点上的死亡人数，然后计算每个死亡点或截尾点上的观察人数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)=poPlP2Pi(11.12)(5) 生存率标准误SEs(ti)的计算，由于s(ti)是依据样本资料计算的，有抽样误差。乘积限估计的标准误的近似计算公式为：rdSEs(ti)=s(ti)(11.13)巧二nj(nj-山)本例s(ti)及标准误的计算列表11.2的(7)、(8)栏中。11.2.2生存率的组间比较在医学随访研究中，通常将病人按随机化方法分配到两种或多种治疗组中，然后随访观察和比较其生存时间的长短和生存率的大小，以此来考察各种治疗方案的优劣；或者分析和比较同一治疗方案下具有不同特征病人的生存率的大小，以此来探讨影响这种疗法的因素。因此，生存率组间比较实际上是两条或多条生存曲线的比较。生存率的假设检验方法有参数法和非参数法两类。参数法要求生存时间已知服从于某种概率分布，对实际资料拟合分布并求得其相应的参数，然后通过比较不同组的分布参数来比较生存率是否相同。非参数法对资料的分布没有要求，适用面比较广。常见的有Logrank检验、Wilcoxon检验(Gehan检验)和似然比检验，似然比检验要求资料服从指数分布才有效。这里主要介绍Logrank检验和Wilcoxon检验两种方法。1. 对数秩和检验(logranktest)Logrank检验是Mantel等人在1966年提出的，这种方法是在组间生存率相同的检验假设(H。)下，对每组生存数据依据在各个时刻尚存活的患者数和实际死亡数计算期望死亡数，然后将期望死亡数与实际死亡数进行比较，作假设检验。这种方法可适合两组或多组生存率比较。这种方法在两组生存率比较时，计算比较简单，故结合实例对两组比较的计算方法作为详细介绍，然后推广到多组检验方法。例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)栏中。计算各组在时间ti上期望死亡数egi。依据H：st1)=S2(t2),egi的计算公式为：(11.14)egi=ngidi/nil求各组实际和期望死亡数的和，并计算X2统计量。各组实际死亡数的和为Dg=Vdgi,i丄l期望死亡数的和为Eg=瓦egi（lX20。05,1=3.84P0.05,两组样本生存率间差异有统计学意义。用猪苓提取物合用化疗治疗急性白血病的生存期大于对照疗法。表11.3对数秩和检验计算表编号i时间用猪苓组（第1组）对照组（第2组）合计ti(1)n1i(2)d1i(3)W1i(4)eii(5)V1i(6)n2i(7)d2i(8)W2i(9)e2i(10)V2i(11)ni(12)di(13)11.516000.6150.236710100.3850.23672612+216010.00.09010.00.025033.515000.6520.22688100.3480.22682314415100.6820.21697000.3180.21692215614000.6670.23817100.3330.23812116+614020.00.06000.00.020076.512000.6670.22226100.3330.22221818+6.512000.00.05010.00.017097.5+12010.00.04000.00.0160108.511100.7330.19564000.2670.195615111+910010.00.04000.00.014012109100.6920.21304000.3080.213013113118000.6670.22224100.3330.22221211411+8000.00.03010.00.01101512+8010.00.30252000.00.010016137101.5560.12242100.4440.30259217176000.8570.01100.1430.12247118186101.0000.00000.06119+195010.00.00000.05020244101.0000.00000.04121263101.0000.00000.03122312101.0000.00000.0212343+1010.00000.010合计811.7882.196473.2122.1964精确法x2统计量的计算公式为：22X=(Dg-Eg)/Vg(11.佝式中Vg为第g组的Fg的方差估计值，其计算公式为：(11.17)ngi5i-ngjdj(nj-djn2(ni-1)所计算的V1=V2，精确法X2值服从自由度为1的X2分布。l本例先求Vgi,再求Vg=v、.gi，得V1=V2=2.1964,x2统计量为：yX2=(811.788)=6.533专业结论同近似法。2.1964对于多组生存率比较，其近似法X2统计量的计算方法与两组相同。设有m组g=1,2,m这时自由度U=m-1。关于多组Logrank检验的精确法x2统计量计算较为复杂，其计算公式为：2_1X=sxVsU=m-1(11.18)其中S，=(S1,S2,s，,m-1),s为向量s的转置。Sg的计算公式为：lSg=(dgi_dingi/nJ=Og-Eg,gh,m-1(11.19)i1V为(m-1)(m-1)矩阵，记为V=Vgh(m-1)(m-1)。其中Vgh为第g组与第k组之间的方差与协方差，计算公式为：lVgh=(nmh、gh-nihnig)di(n-dJ/nf(n-1)(11.20)i=1当g=h时，Sgh=1;当gzh时，Sgh=0。从多组生存率比较的计算公式可知，当m=2时，即为两组生存率的计算公式。Wilcoxon检验法及与Logrank检验法比较Wilcoxon检验法当g=1,2,g,伸寸,Wilcoxon检验法x2统计量计算公式仍可表示为：x=svsu=m-1其中s=(S1,S2,s，,m-1),S，为向量s的转置。Sg的计算公式为：lsg=Wi(dig-dinig/nJ(11.21)i巧V为(m-1)(m-1)矩阵，记为V=Vgh(m-1)(m-1)Vgh的计算公式为：l22Vgh=二wi(ninih：；gh-nihnig)di(ni-di)/ni(ni-1)(11.22)i4上面Sg和Vgh计算公式中Wi为权重，这里Wi=ni。对于例11.3资料，若用Wilcoxon检验法，得x2=5.2822，结论同Logrank检验。(2)两种检验方法的比较Logrank检验法和Wilcoxon检验法实际上可以用统一公式来表示，即Wilcoxon检验法的公式，公式中的权重Wi=1时为Logrank检验法，wi=ni时为Wilcoxon检验。因而可以发现，Logrank检验对生存时间较长的个体在检验中权重较大，对生存时间较短的个体在检验中权重较小，在生存率(曲线)比较中，这种方法对尾部较为敏感；而Wilcoxon检验则与Logrank检验相反，对生存时间较短的个体在检验中权重较大，比较中对数据的头部的差别较为敏感。从理论和实践中均发现，当生存资料的各死亡点的危险率在两组或多组间成比例时，Logrank检验的效率高于Wilcoxon检验法，宜选用Logrank检验。当生存资料各时点的危险率服从其他状态时，Wilcoxon检验法效率高于Logrank检验，宜选用Wilcoxon检验法。11.2.3应用实例例11.4对例11.3的资料，试估计两组的生存函数，并用Logrank检验和Wilcoxon检验法作比较。1. SAS程序：DATASA;INPUTT;IFT16THENGROUP=2;ELSEGROUP=1;T=ABS(T);CARDS;-24-6-6-7.58.5-910-121318-19242631-431.5-23.566.5-6.511-111317PROCLIFETESTMETHOD=PLPLOT=(S);TIMET*CENSOR(1);STRATAGROUP;RUN;程序说明数据输入时截尾数据用负值表示，第一个if语句产生指示变量censor，其取值为1时为截尾数据，取值为0时为完全数据。第2个if语句产生分组变量group，前16个数值属于用猪苓组(1组)，后10个数据属于对照组(2组)。对t取绝对值是为了保证参与计算的生存时间都是正值。PROCLIFETEST过程可选择寿命表法或乘积限估计法对生存资料作生存率估计，并可对两组或多组间的生存率作Logrank检验和Wilcoxon检验。method=选项有两种选择，即PL和life，PL选项为隐含值，表示用乘积限估计法估计生存率，Life表示用寿命表法估计生存率，当用寿命表法分析时，程序自动形成生存时间的分组区间，也可人为指定分组区间，在proc语句中加上intervals=(atobbyc),a、b、c分别表示初值、终值和步长。Plots=()要求绘图，其中s表示生存率曲线，L表示取对数，H表示危险率函数，图形的横纵坐标分别为：s=(t,s)、LS=(t,-log(s),LLs=(log(t),log,Log(-log(s)、H(t,h)time语为lifetest过程中必须语句，为设置生存时间和截尾指示。当有分组变量时用strate表示。如果资料中还包含有协变量，可在num语句前加test语句，如testX1X2，以便检验协变量与生存时间联系的密切程度。如果用下一节介绍的phreg过程来揭示协变量与生存时间的关系则更好。输出结果及解释结果首先输出两组的乘积限法估计的生存率(survival)及死亡率(Failure)，生存率的标准误(SurvivalStandardError)、死亡数(NumberFailed)和存活数(Numberleft)，表示有*号者为截尾观察值。输出中有生存时间的四分位数，第50%位数为中位生存期，第1组为24个月，第2组为11个月，两组平均生存期分别为20.35和10.76个月，由于第1组最后一个值为截尾数据，所以均数的估计是有偏性的。TheLIFETESTProcedureProduct-LimitSurvivalEstimatesGROUP=1SurvivalStandardNumberNumberTSurvivalFailureErrorFailedLeft0.00001.0000000162.0000*0154.00000.93330.06670.06441146.0000*1136.0000*1127.5000*1118.50000.84850.15150.09992109.0000*2910.00000.75420.24580.12563812.0000*3713.00000.64650.35350.14684618.00000.53870.46130.15705519.0000*5424.00000.40400.59600.16576326.00000.26940.73060.15597231.00000.13470.86530.12318143.0000*80*CensoredObservationSummaryStatisticsforTimeVariableTQuantileEstimateLower,Upper)75%31.000018.000050%24.000013.000031.000025%13.00008.500024.0000Mean20.3549StandardError2.9640NOTE:Thelastobservationwascensoredsotheestimateofthemeanisbiased.TheLIFETESTProcedureProduct-LimitSurvivalEstimatesGROUP=2SurvivalStandardNumberNumberTSurvivalFailureErrorFailedLeft0.00001.0000000101.50000.90000.10000.0949192.0000*.183.5000*.176.00000.77140.22860.1442266.50000.64290.35710.1679356.5000*.3411.00000.48210.51790.18774311.0000*.4213.00000.24110.75890.19465117.000001.0000060*CensoredObservationCensoredObservationSummaryStatisticsforTimeVariableTPoint95%ConfidenceIntervalQuantileEstimateLower,Upper)75%13.000011.000017.000050%11.00006.000017.000025%6.50001.500013.0000Mean10.7571StandardError1.9434SummaryoftheNumberofCensoredandUncensoredValuesGROUPTotalFailedCensored%Censored12Total168850.0000106440.000026141246.1538以下输出Kaplan-Meier生存曲线图TheLIFETESTProcedureSurvivalFunctionEstimatesSDF|1.0+*-*-A|AAS|BB|ru|A|-Aiv0.8|+|v|BB|laD|B-A-|-|-B-AAAD|BBAAis0.|6+|tr|AAib|B-B|i0.4+o|nF|F|u|n|c0.2+itco0.|2|+o|n|0.0+|B|BBB|A-A-|A-+051015202530354045TStrataB+BBBBA+A*AAAAA-+051015202530354045TLegendforStrataSymbolsA:GROUP=1B:GROUP=2最后输出Logrank检验和Wilcoxon检验过程中的秩统计量(x2统计量和p值，其结果与例11.3的手工计算结果一致。TestingHomogeneityofSurvivalCurvesoverStrataTimeVariableTRankStatisticsGROUPLog-RankWilcoxon-3.1355-46.0003.135546.000CovarianceMatrixfortheLog-RankStatisticsGROUP121.95379-1.95379-1.953791.95379CovarianceMatrixfortheWilcoxonStatisticsGROUP121577.500-577.500S）和方差协方差（V）以及2-577.500577.500TestofEqualityoverStrataPrTestChi-SquareDFChi-SquareLog-Rank5.032110.0249Wilcoxon3.664110.0556-2Log(LR)2.223910.135911.3Cox回归模型与PHREG过程1972年,英国统计学家DR.Cox提出了半参数生存分析数学模型Cox回归模型（Coxregressionmedel）,在以后二十多年中，众多的生物统计学家在理论和应用方面作了大量的研究。目前，Cox回归模型已成为生存分析中理论较为完善，应用最广泛的统计模型。我国随着计算机和统计软件的推广应用也已开始应用Cox模型进行慢性病的预后分析。实践证明，Cox模型对许多生存资料都有用，而且行之有效。本节简要介绍Cox模型结构、意义、参数估计和假设检验方法，以及SAS软件中如何应用PHREG过程作Cox模型分析。11.3.1Cox模型的结构及流行病学意义在随访研究中，要考察和比较不同的治疗方法，不同病理类型，病人的某些特征对疾病预后的影响，可运用非参数方法进行组间比较，这类方法虽然使用简单，但在多因素共存条件下，单一因素的比较会受到其他因素干扰和混杂，组间难以达到均衡，而且不能分析和考察因素间的关系(如交互作用)和进行定量评价。因此，需用多因素回归模型的方法。生存分析中，最典型、最常用的就是Cox回归模型。Cox模型的结构对于一组带有多个预后因素(称为解释变量或协变量)的随访病人，假设病人不受到任何协变量的作用，那末，病人在整个随访期间会显现一定的生存或死亡规律。这一死亡规律用h(t)来描述和表示，而且所有病人的死亡规律应该是一致的。而病人在随访期间实际上是受到多个预后因素的综合作用。预后因素对病人的生存或死亡规律的影响可以看作是协变量对ho(t)的修改。并且病人各自所带有的预后因素是不同的。因而，病人实际上呈现出不同的生存或死亡规律。病人的实际死亡规律用h(t,x)表示(x表示协变量)。由此可见，生存分析回归模型由基本部分ho(t)和修改部分$(x)($(x)为预后因素x的函数)两部分组成。这里协变量X=Xi,X2,，x可以是不同的治疗方法或手术方法，不同病理类型或病人的性别、年龄等特征，也可以是复合变量(如交互作用)。当基本部分ho(t)与修改部分(x)为乘积关系时，即危险函数由基本危险率函数乘上一个协变量的常数因子组成。h(t,x)=ho(t)(x)(11.23)称此模型为成比例危险模型(Proportionalhazardsmodel)。带协变量的修改部分(x)最常见的形式为exp(x3则模型为：h(t,x)=h0(t)exp(x3=ho(t)exp(31X1+32x2+3pxp)(11.24)其中3=(31,32,3p)是p个未知待估计的回归系数，它是描述各个因素对生存期影响大小的参数向量；X=(X1,X2,寿)是一个p维协变量向量，xp可以是影响预后的治疗方法、病人的某些特征，也可以是这些因素的交互作用项；ho(t)是未知的、任意分布的基准危险函数。对比例危险模型，若假定h(t)为特定的时间分布函数，则称为全参数生存分析模型，可以用极大似然法作参数3估计和检验。当h(t)未知分布时，用部分似然函数(Partiallikelihoodfunction)对参数3作出估计，则称此模型为半参数模型，也就是通常意义下的Cox回归模型。流行病学意义应用部分似然函数理论对模型中的3作出估计时，其意义在于：当某个协变量Xj为二分变量时，Xj=1的危险性相对于Xj=0为：(11.25)hoWexpX!+岛.1+BpXp)RRR=hoWexgx!汕pXp)=e从上式可以发现Xj=1相对xj=o在任意时间的危险性是一个常数e，也即两组病人的相对危险性为一个常数，成比例危险模型由此而得名。当为为多分类有序变量或连续变量时，elj为变量每增加一个单位，个体增加或减少的相对危险性，变量Xj是危险因素还是保护因素要视3j的正负值和变量的取值而定。11.3.2Cox模型的参数估计与统计推断1Cox模型的似然函数当生存数据中没有重合死亡(失效)点(无ties)时，所有样本个体(i=1,2,n的生存时间(i=1,2，,d和截尾时间可以看作是同一始点开始随访观察。按从小到大排序记为ti,t2,n,可以看作这一群个体须经过不同的死亡点(Failure)而在各死亡点必须有一个个体死亡(失效)。经过死亡点的所有个体记作死亡点ti的危险集R(ti),已经生存到ti时刻的个体(R(ti)，在ti时刻必须有一个个体死亡，而这个个体就是第i个体的概率，称为条件死亡概率或条件危险率，记为：P(i/R(ti),3)给定R(ti)且已知ti时刻有一个病人死亡的条件下，则病人是第i个体的概率是：(11.26)P(i/R(ti),3)=hi(t,x)/vhk(t,x)心)根据似然函数构造理论，Cox模型的似然函数应是各个死亡点对它的贡献而构成，即各死亡点ti的条件危险率的乘积：dL(3)=1Ii=1hi(t,x)hk(t,x)=丨丨i生ho(t)i(x)hk(t);：k(x)d=11idkR(tJkR(tl)exp(xi：)、exp(Xk：)kR(tl)有趣的是在有截尾数据时，其似然函数仅为死亡病人在各死亡点上条件死亡概率的乘积。这一似然函数并不是通常意义下的似然函数。DRCox于1975年证明了这是特定意义下全似然函数中的一部分，故称部分似然函数(Partiallikelihoodfanction)。2.参数3及其方差的估计Cox也证明了可将一般极大似然理论应用到部分似然函数的参数3及其方差的估计。为计算方便，对Cox模型的部分似然函数两边取对数，得对数部分似然函数。d(11.27)lnL(3)=二xi：-1n二exp(xk：)i=1kR(ti)应用极大似然估计理论，将对数部分似然函数即：lnL(3)对3j求一阶偏导数并令其等于0。汕L(-)_n1=0TnL(-)=0P个非线性方程组；-2a:TnLJ门=0帝0对这一似然方程组，可用迭代法(例如Newton-Raphson法)求3j的估计值(迭代求解过程从略)。参数3j的估计值?j的方差var(?j)和标准误SE(?j)可用二阶偏导数I讥)/川工1一的负值组成的信息矩阵|pp及其逆矩阵l；p二Vpp(方差协方差矩阵)的主对角线元素求得。3参数3及回归效果检验在cox模型的参数估计和模型建立过程中，要对变量进行筛选并作出检验。DRCox的1972和1975年的文章以及以后十多年中许多学者对部分似然函数的大样本性质作了证明，提出部分似然可代替普通似然对参数3进行统计推断。关于极大部分似然估计3的统计推断有三种渐近方法。似然比检验(Likelihoodratiotest)设原模型参数为3=(31,32,p),并求得对数部分的似然函数值InL(P)；需检验的参数3=(31,32,，3k)(kPROCPHREG;MODELT*D(0)=GROUPRENAL/TIES=BRESLOWSELECTION=SW;RUN;2.SAS程序说明建立数据集时第一列为生存天数，第二列为截尾指示，截尾数据为0，生存数据为1。第三列为组别，A组为1,B组为2。第四列为肾功能损害情况，有损害为1,未损害为0。程序调用PHRE（过程分析，PHREG为proportionalhazardregression（比例风险回归）。PHREGMODE为Cox模型分析的必须语句，MODE语句左边为生存时间和截尾指示，右边为需分析的解释变量。选择项ties=Breslow|Discrete|Efron|Exact，指定重复失效时间数据的处理方法。Breslow为隐含值，当有少量重合数据时，程序自动选用1974Breslon提出的近似似然；Discrete为选用离散Logistic模型，在重合数据较多或生存时间为离散时可选此项；Exact为Breslow和Efron的精确算法。当需对解释变量进行筛选时，应选用selection=FW|BW|SW|score和最优子集选择法选项，分别表示向前法、向后法、逐步法和最优子集选择法。引入或剔除变量的显著性水平可通过sle=和sls=来指定，隐含值a=0.15。若要对资料进行分层cox模型分析，可在run语句前加STRATA语句指定分层变量。输出结果及解释ThePHREGProcedureDataSet:WORK.COXDependentVariable:TCensoringVariable:DCensoringValue(s):0TiesHandling:BRESLOWSummaryoftheNumberofSummaryoftheNumberofEventandCensoredValuesPercentTotalEventCensoredCensored2521416.00Step1:VariableRENALisentered.Themodelcontainsthe