事件史分析操作及说明A＿08

事件史分析的Stata操作及说明 目录 离散时间模型示例：职位晋升分析...................................................................................................................... 1 Cox模型 ................................................................................................................................................................... 2 Cox模型示例：癌症数据分析 ............................................................................................................................... 3 分层的(stratified) Cox模型的背景 ....................................................................................................................... 11 示例：癌症生存的分层分析................................................................................................................................ 12 最大似然估计的连续时间存活分析.................................................................................................................... 18 Exponential分布模型 ............................................................................................................................................ 18 Gompertz模型........................................................................................................................................................ 19 Weibull模型 ........................................................................................................................................................... 19 用Stata拟合exponential模型、Weibull模型和Gompertz模型：癌症生存 ........................................................ 19 本章参考文献： 1. * Lawrence C. Hamilton. Statistics with STATA, Updated for Version 9. Brooks/Cole, a division of Thomson Learning, Inc.,-2006. （第11章）// 中译本：郭志刚等，应用STATA做统计分析，重庆大学出版社，2008. 2. * Rabe-Hesketh and Everitt. A Handbook of Statistical Analyses Using Stata. Chapman & Hall/CRC, 2004 3. * Kleinbaum. Survival Analysis. 1996; 4. * Yamaguchi, Kazuo. Event History Analysis. Sage Publications. Inc. 1991 5. * Allison, Paul D.. Event History Analysis: Regression for Longitudinal Event Data. Sage Publications, Inc. 1984 6. * Blossfeld, Hamerle and Mayer. Event History Analysis. Lawrence Erlbaum Associates, Inc. 1989 7. * Blossfeld and Rohwer. Techniques of Event History Modeling: New Approaches to Causal Analysis (Second Edition). Lawrence Erlbaum Associates, Publishers, 2002. 8. * Cleves, Gould and Gutierrez. An Introduction to Survival Analysis Using Stata (reviced edition). Stata Press Publication, 2004. 9. * Janet M. Box-Steffensmeier and Bradford S. Jones. Event History Modeling : A Guide for Social Scientists. Cambridge University Press, 2004 10. 梁在，事件史分析，载郭志刚主编《社会统计分析方法----SPSS软件应用》，中国人民大学出版社，1999年。 11. 郭申阳，使用SPSS软件对事件史原始数据进行预处理，载郭志刚主编《社会统计分析方法----SPSS软件应用》， 中国人民大学出版社，1999年。 12. 郭志刚，历时研究与事件史分析，《中国人口科学》第1期，2001年2月，第67～72页。 13. Cox, "Regression models and life tables," Journal of the Royal Statistical Society, 1972 14. Petersen, "The statistical analysis of event histories," Sociological Methods and Research, 1991 15. Allison, "Discrete-time methods for the analysis of event histories," Sociological Methodology, 1982 16. Poston, "Son preference and fertility in China." Journal of Biosocial Science 34: (2002): 333-347 标志 * 者，社会学系资料室有正版或复印本。 重要声明 本讲义仅作为选修《高级社会统计专题》课同学的参考资料。 本讲义内容是从有关文献中选取、翻译和编辑的，辅助学生理解统 计方法原理和操作。本讲义不允许复印、传播或用于其他目的。 使用中如发现错误之处，请联系：zguo@ pku.edu.cn D:\zguo_x60\统计及SPSS课程\社会统计课程立项\上网\6.生存分析_1\事件史分析操作及说明A＿08.doc 2008-10-27 郭志刚·北京大学社会学系 1 事件史分析的Stata操作及说明 离散时间模型示例：职位晋升分析 用梁在（1999：p400，例 1）关于职位晋升的数据来示范。该数据原为 SPSS 系统数据文件，已 经用 StatTransfer 软件转换为 Stata 系统数据文件了，文件名为 T12_3.dta。该数据有 5 个变量：id（案 例识别码）、dura（在本公司工作时间）、prom（该人年中是否晋升）、wexp（到本公司前是否有工 作经历）、sex（性别）。用 dura、wexp 和 sex 来解释 prom。 应用 Stata 中命令形式为： logit Y x1 x2 x3 其中 Y 代表虚拟因变量，各 X 代表自变量。对于本例的具体命令为： logit prom dura wexp sex Logit estimates Number of obs = 52 LR chi2(3) = 10.85 Prob > chi2 = 0.0126 Log likelihood = -28.712942 Pseudo R2 = 0.1589 ------------------------------------------------------------------------------ prom | Coef. Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- dura | 1.5433 .7046274 2.19 0.029 .1622553 2.924344 wexp | 1.726036 .7677605 2.25 0.025 .2212527 3.230819 sex | 1.2279 .6775754 1.81 0.070 -.1001238 2.555923 _cons | -4.00691 1.358362 -2.95 0.003 -6.669249 -1.34457 ------------------------------------------------------------------------------ 上述输出与 SPSS 输出结果相同。回归结果表明，在本公司工作时间 dura 和以往工作经历 wexp 都对晋升有显著的正影响。虽然性别 sex 的回归系数为正值，表现出男职员比女职员更有可能晋升， 然而该变量并未显著。 在 logistic 回归时，我们更愿意用优势比或发生比率(Odds Ratio)的方式来解释自变量的影响作 用。优势比＝e(b) 。用下列 Stata 命令可以直接输出优势比。 logistic prom dura wexp sex Logistic regression Number of obs = 52 LR chi2(3) = 10.85 Prob > chi2 = 0.0126 Log likelihood = -28.712942 Pseudo R2 = 0.1589 ------------------------------------------------------------------------------ prom | Odds Ratio Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- dura | 4.680007 3.297661 2.19 0.029 1.17616 18.622 wexp | 5.618337 4.313537 2.25 0.025 1.247639 25.30036 sex | 3.414051 2.313277 1.81 0.070 .9047254 12.88319 ------------------------------------------------------------------------------ 优势比系数表明，在其他条件不变时，在本公司工作时间每增加 1 年，晋升发生比将为原来的 4.7 倍。进入本公司前有工作经历的人的晋升发生比是以前没有工作经历的人的 5.6 倍。男职员的晋 D:\zguo_x60\统计及SPSS课程\社会统计课程立项\上网\6.生存分析_1\事件史分析操作及说明A＿08.doc 2008-10-27 郭志刚·北京大学社会学系 2 升发生比是女职员的 3.4 倍，但统计性没有达到显著 标准 excel标准偏差excel标准偏差函数exl标准差函数国标检验抽样标准表免费下载红头文件格式标准下载 。 也就是说，实际上在用 logistic 回归分析事件风险时，其自变量影响的解释与常规 logistic 回归一 样。所不同的只是现在用的是可以表现历时风险变化的人年 记录 混凝土 养护记录下载土方回填监理旁站记录免费下载集备记录下载集备记录下载集备记录下载 数据。 Cox模型 Cox 模型是人口学和社会学中应用最普遍的一种事件史分析模型，全称为 Cox 比例风险模型 （Cox proportional hazards method.）。 Yamaguchi 在其事件史著作中甚至都没有讨论最大似然估计的连续时间模型。他直接使用 Cox 提 出的偏似然估计方法，并将此称为连续时间的比例风险模型。他还认为，Cox 模型是应用最普遍的事 件史分析模型（p101）。 第一，Cox 模型最显著的优点是不需要定义具体形式就可以假定生存时间的依赖关系。在后面介 绍那三种 ML 连续时间模型时就可以看到，研究者必须在模型中定义具体的时间分布形式。但是在 Cox 模型中，研究者并不需要定义具体时间分布形式。 第二，Cox 模型具有采用分层(stratified)分析的能力。这里，分层分析指模型可以控制某一分类 协变量中的某一类、或某一组分类协变量的影响。它们的影响可能与时间之间存在复杂的互动关系。 Cox 比例风险模型是 ML 模型的推广。当模型中有两个固定自变量、且没有动态自变量时，Cox 模型可以用下列 公式 小学单位换算公式大全免费下载公式下载行测公式大全下载excel公式下载逻辑回归公式下载 表达： h t = a t b X b X1 1 2 2log ( ) ( )+ + 其中a t( )是时间的任意函数。但是，这一时间函数并没有具体定义。这就是 Cox模型的诱人之 处。因此，Cox模型又被称为偏参数模型或半参数模型（partially parametric 或 semi-parametric）， 因此时间函数不一定要被定义。 Cox 模型是用偏似然估计求解的，这使其具有天然优势，因为偏似然估计既能容许时间函数存 在，又不需要具体定义具体函数形式。 比例风险模型假定，风险率为协变量影响参数的对数线性函数。用 ( )ih t 表示个人 i在时点 t的 风险率，其公式为： i 0 k i kh t = h t b X t ( ) ( ) exp[ ( ) ]× (式1) 简称上式为式 1。其中： h t0 ( ) 代表时间函数的主要特征，称为基准风险函数； i kX t( ) 是个人 i在时点 t上的第 k个自变量，既可以依时间变化，也可以是固定的。 基准风险函数h t0 ( )为所有案例所共有。但是，在这一方程用 Cox 的 PL 估计方法求解时，h t0 ( ) 的函数形式并不用具体定义。我们来看看这是怎么回事（参见 Yamaguchi：p106-107 和 Allison：p69- 71）。 参数的 PL 估计通过使偏似然函数最大化来求出。PL 函数的建立过程为：首先，在持续期 (duration，即风险期)长度的基础上，将案例按其持续期由小到大排序。下面公式中的下标 i将代表 重新排序后的案例编号。于是，建立的 PL函数为： i I i i j i i j i PL h t h t 1 [ ( ) / ( )]δ = ≥ = ∏ ∑ (式2) 其中： j ih t( )为第 j个案例在时点 ti上的风险函数值，其中 ti是第 i个案例发生事件或被删截 的时点。式 2 中的δi是一个虚拟变量，当第 i个案例发生事件时为 1，而发生删截时为 0。Π是表示连 D:\zguo_x60\统计及SPSS课程\社会统计课程立项\上网\6.生存分析_1\事件史分析操作及说明A＿08.doc 2008-10-27 郭志刚·北京大学社会学系 3 乘的符号。 当把式 1中的 ih t( )代入式 2 的分子和分母中时，其中的基准风险函数h t0 ( )就被约分掉了。于 是，PL 函数可以仅仅表达为协变量参数的函数，如下： i I ik jkk k ik ki i j i PL tb b tX X 1 {exp[ )]/ exp[ ( )]}( δ = ≥ = ∏ ∑ ∑ ∑ (式3) 上述过程意味着，尽管在 Cox模型中存在基准风险函数，它反映出风险率的时间分布形式（见 式 1），然而当模型实际中是用 PL方法来估计的时候，其具体函数形式并不需要定义，因为不管是 什么函数形式，都会同时被从式 2的分子和分母中被删除。所以，最终方程（式 3）中并不存在基准 风险函数h t0 ( )。 正是因为 Cox 模型不需要具体定义风险的时间分布形式，所以它的应用要容易得多，也应用得 很普遍。 Cox模型示例：癌症数据分析 下面用 Stata 软件附带的 cancer_01.dta 数据来示范估计 Cox 模型。我们先来熟悉一下这个数据。 use "cancer_01.dta", clear (Patient Survival in Drug Trial) desc Contains data from C:\Stata\cancer.dta obs: 48 Patient Survival in Drug Trial vars: 4 size: 576 (99.9% of memory free) ------------------------------------------------------------------------------- studytim int %8.0g Months to death or end of exp. died int %8.0g 1 if patient died drug int %8.0g Drug type (1=placebo) age int %8.0g Patient's age at start of exp. drug01 float %9.0g 1=real drugs, 0=placebo agemean float %9.0g age centered at its mean ------------------------------------------------------------------------------- 数据中原来仅有前 4 个变量。表达生存的因变量实际上包括其中的前两个，后两个变量将作为协 变量。第 5 个变量是将第 3 个变量改造为一个简单的虚拟变量。因为少设了一个虚拟变量（旨在示范 时简单），改造过程损失了一些信息。 各变量的意义如下： 1、 studytim 是癌症患者的风险期长度，即自患病始至死亡或删截的时间。观测值域为 1~39 个月。 2、 died 是个虚拟变量，1 表示风险期以死亡结束（即无删截，占 65%），0 表示删截。 3、 drug 是原始分类变量，表示患者采用的药物，1 是安慰剂，2 和 3 代表两种真实的不同药物。 4、 age 是患者的年龄，按岁计量。 5、 drug01 是使用药物的虚拟变量：0 是安慰剂，1 为使用真实的不同药物（未区分具体药物）。 6、 agemean 是将 age 以其平均值对中的改造变量。（该变量在做分层的生存曲线图时有用。） 在用 Stata 做 EVA 时，必须先要用 stset 命令（表示 survival time set）声明将要用这个数据做生 存分析。随后要依次定义生存时间的变量 studytim 和表示这一风险期是否删截的变量 died。（注意这 D:\zguo_x60\统计及SPSS课程\社会统计课程立项\上网\6.生存分析_1\事件史分析操作及说明A＿08.doc 2008-10-27 郭志刚·北京大学社会学系 4 样定义以后 Stata 将自动转换数据格式以适应生存分析的需要！） stset studytim died failure event: died != 0 & died < . obs. time interval: (0, studytime] exit on or before: failure ------------------------------------------------------------------------------ 48 total obs. 0 exclusions ------------------------------------------------------------------------------ 48 obs. remaining, representing 31 failures in single record/single failure data 744 total analysis time at risk, at risk from t = 0 earliest observed entry t = 0 last observed exit t = 39 下面要求输出生存数据的概要统计。 stsum failure _d: died analysis time _t: studytime | incidence no. of |------ Survival time -----| | time at risk rate subjects 25% 50% 75% ---------+--------------------------------------------------------------------- total | 744 .0416667 48 8 17 33 输出表明，有 48 个患者，共计 744 个月生存期。如果假定每个月的死亡风险率都相同的话，那 么可以估计出平均每个人月的死亡风险率为 0.0417，即死亡 31 人与 744 个人月之比值。 上表中failures依题意就表示死亡，744个人月中发生死亡的频数为31，因此stsum表中的死亡风险为 31/744 = .04167（这其实就是所谓的风险函数值，即发生率incidence rate）。生存时间的四分位值是 从K-M（Kaplan-Meier）生存函数（见后面）中推导出来。从K-M生存函数估计出，死于前8个月之内 的平均概率为25%，死于前17个月内的平均死亡概率为50%（提示：不要将这17个月作为中位存活时 间和死亡时间来理解！1），死于前33个月内的死亡概率达到75%。 下面取得这一生存数据的描述统计。 stdes failure _d: died analysis time _t: studytime |-------------- per subject --------------| Category total mean min median max ------------------------------------------------------------------------------ no. of subjects 48 no. of records 48 1 1 1 1 (first) entry time 0 0 0 0 (final) exit time 15.5 1 12.5 39 subjects with gap 0 time on gap if gap 0 time at risk 744 15.5 1 12.5 39 1 它是年龄别生存概率的中位数，而不是年龄别尚存人数的中位数。 D:\zguo_x60\统计及SPSS课程\社会统计课程立项\上网\6.生存分析_1\事件史分析操作及说明A＿08.doc 2008-10-27 郭志刚·北京大学社会学系 5 failures 31 .6458333 0 1 1 ------------------------------------------------------------------------------ 输出表明，有 48 个患者，共计 744 个月生存期，其中已经死亡 31 人。 从 K-M 生存函数估计出，患者平均生存 15.5 个月。最短的只生存了 1 个月，中位生存月数为 12.5 个月。（注意：此值不等于前表输出中的平均死亡概率达到 50%时所对应的时间！）最长的生 存期达到 39 个月。已经死亡的 31 个人占总数 48 人的 64.58%。 用常规方法取得概要统计，即 sum Variable | Obs Mean Std. Dev. Min Max -------------+-------------------------------------------------------- studytime | 48 15.5 10.25629 1 39 died | 48 .6458333 .4833211 0 1 drug | 48 1.875 .8410986 1 3 age | 48 55.875 5.659205 47 67 _st | 48 1 0 1 1 -------------+-------------------------------------------------------- _d | 48 .6458333 .4833211 0 1 _t | 48 15.5 10.25629 1 39 _t0 | 48 0 0 0 0 这一统计输出证明了对生存期的描述统计的理解是对的，即 48 个人平均生存 15.5 个月，最小值 为 1 个月，最大值 39 个月。死亡人数占到约 65%，换句话说就是有 35%的案例是以删截结束。注意 drug 为 1 时表示安慰剂，2 和 3 才代表真实药物。那么这个变量的平均数 1.875 其实没有意义。所 以，下面我们需要专门对 drug 列出频数表来看看分布。变量 age 的均值为 55.9 岁，最小者 47 岁，最 大者 67 岁。我们还看到 Stata 在将数据转换为生存时间格式时还新加了 4 个变量。 下面我们要求列出 drug 的频数表： tab drug Drug type | (1=placebo) | Freq. Percent Cum. ------------+----------------------------------- 1 | 20 41.67 41.67 2 | 14 29.17 70.83 3 | 14 29.17 100.00 ------------+----------------------------------- Total | 48 100.00 变量 drug 的频数表说明，有 42%的案例使用安慰剂，29%用的是第 1 种药物，另外有 29%用的 是第 2 种药物。 tab drug drug01 Drug type | (1=placebo | drug01 ) | 0 1 | Total -----------+----------------------+---------- 1 | 20 0 | 20 2 | 0 14 | 14 3 | 0 14 | 14 D:\zguo_x60\统计及SPSS课程\社会统计课程立项\上网\6.生存分析_1\事件史分析操作及说明A＿08.doc 2008-10-27 郭志刚·北京大学社会学系 6 -----------+----------------------+---------- Total | 20 28 | 48 以上可看出 drug 是如何虚拟变量化的。 由于数据仅有 48 个案例，不妨用常规方法列出所有案例的变量值。 list +----------------------------------------------------+ | studyt~e died drug age _st _d _t _t0 | |----------------------------------------------------| 1. | 1 1 1 61 1 1 1 0 | 2. | 1 1 1 65 1 1 1 0 | 3. | 2 1 1 59 1 1 2 0 | 4. | 3 1 1 52 1 1 3 0 | 5. | 4 1 1 56 1 1 4 0 | |----------------------------------------------------| 6. | 4 1 1 67 1 1 4 0 | 7. | 5 1 1 63 1 1 5 0 | 8. | 5 1 1 58 1 1 5 0 | 9. | 8 1 1 56 1 1 8 0 | 10. | 8 0 1 58 1 0 8 0 | |----------------------------------------------------| 11. | 8 1 1 52 1 1 8 0 | 12. | 8 1 1 49 1 1 8 0 | 13. | 11 1 1 50 1 1 11 0 | 14. | 11 1 1 55 1 1 11 0 | 15. | 12 1 1 49 1 1 12 0 | |----------------------------------------------------| 16. | 12 1 1 62 1 1 12 0 | 17. | 15 1 1 51 1 1 15 0 | 18. | 17 1 1 49 1 1 17 0 | 19. | 22 1 1 57 1 1 22 0 | 20. | 23 1 1 52 1 1 23 0 | |----------------------------------------------------| 21. | 6 1 2 67 1 1 6 0 | 22. | 6 0 2 65 1 0 6 0 | 23. | 7 1 2 58 1 1 7 0 | 24. | 9 0 2 56 1 0 9 0 | 25. | 10 0 2 49 1 0 10 0 | |----------------------------------------------------| 26. | 11 0 2 61 1 0 11 0 | 27. | 13 1 2 62 1 1 13 0 | 28. | 15 0 2 50 1 0 15 0 | 29. | 16 1 2 67 1 1 16 0 | 30. | 19 0 2 50 1 0 19 0 | |----------------------------------------------------| 31. | 20 0 2 55 1 0 20 0 | 32. | 22 1 2 58 1 1 22 0 | 33. | 23 1 2 47 1 1 23 0 | 34. | 32 0 2 52 1 0 32 0 | 35. | 6 1 3 55 1 1 6 0 | |----------------------------------------------------| 36. | 10 1 3 54 1 1 10 0 | 37. | 17 0 3 60 1 0 17 0 | 38. | 19 0 3 49 1 0 19 0 | 39. | 24 1 3 58 1 1 24 0 | 40. | 25 0 3 50 1 0 25 0 | |----------------------------------------------------| D:\zguo_x60\统计及SPSS课程\社会统计课程立项\上网\6.生存分析_1\事件史分析操作及说明A＿08.doc 2008-10-27 郭志刚·北京大学社会学系 7 41. | 25 1 3 55 1 1 25 0 | 42. | 28 1 3 57 1 1 28 0 | 43. | 28 0 3 48 1 0 28 0 | 44. | 32 0 3 56 1 0 32 0 | 45. | 33 1 3 60 1 1 33 0 | |----------------------------------------------------| 46. | 34 0 3 62 1 0 34 0 | 47. | 35 0 3 48 1 0 35 0 | 48. | 39 0 3 52 1 0 39 0 | +----------------------------------------------------+ 注意第 1 个和第 2 个患者都仅仅参加研究 1 个月，并且也都没有删截，即他们都已经死亡了。第 10 个患者参加研究 8 个月，但是以删截结束。第 45 个患者参加研究 33 个月后死亡了，第 46、47、 48 个患者生存期分别为 34、35、39 个月，并且都仍然生存，即他们都以删截结束。 我们也看到了 stset 新加的 4 个变量。_st 在此例中取值都是 1，是个常量，其意义可能是时间 变量值都不缺失，结果导致_t 值就等于 studytime 值。并且_d 值就等于 died 值，那么意义就是标志 时间变量是否删截。_t0 都等于 0，这标志时间的起点。 我们知道这些患者参加研究的起始时间可能实际上不同，观测时间到可能是一样的。但是现在将 始点都定为 0，然后用观测时间作为时间刻度。这好比： 观测时点 ↓ Î ↓ _t0 SSSSSSSSSSX SSSSSSSSSSX SSSQ SSSQ SSSSSSSSSSSSSSX SSSSSSSSSSSSSSX 下面介绍所谓的 Kaplan-Meier 生存函数。 令 nt代表尚未失败（本例即死亡）且在时期 t开始时未删截的观测案例数，dt代表在时间 t中 的失败案例数（本例即死亡数）。Kaplan-Meier 对超出时期 t的生存人数估计为时期 t的生存概率 与生存时间（人年数或人月数）的乘积： t j j j j t S(t) = n d n 0 {( ) / } = −∏ mdn mdn mdn x x x x mdn mdn 0 mdn s q p l p l =1 l 0 0 0 1 (1 ) = = − = − = = × = ∏ ∏因为： 其中： x mdn x =x when l l0( ) ( 0.5 )=所以： 对上述癌症数据而言，有 2 名患者在 1 个月后就死亡了。在这么短的时间里尚未发生删截，所以 超出时期 t=1的生存概率为： S(1) = (48-2) / 48 = 0.9583 即：(1- q1) 此式就是存活概率（1-qx）：分母为期初人数，分子为期末时（即下期初）的存活人数。 第 1 个时期结束时的生存概率等于本期存活概率。 在时期 t=2中，有 1 人死亡。在时期 t=3中，又有 1 人死亡。因此： S(2) = 0.9583 × (46-1) / 46 = 0.9375 即：S(1)(1-q2)=(1-q1) (1-q2) S(3) = 0.9375 × (45-1) / 45 = 0.9167 即：S(2)(1-q3)=(1-q1) (1-q2) (1-q3) 人口统计学原理告诉我们，当这个生存函数的初始值为实际最初人数时，生存函数值等于 t 期末 时的存活人数（即生命表中的 lx，当 l0＝1 时）。当这个生存函数的初始值为 1.0 时，生存函数值等 D:\zguo_x60\统计及SPSS课程\社会统计课程立项\上网\6.生存分析_1\事件史分析操作及说明A＿08.doc 2008-10-27 郭志刚·北京大学社会学系 8 于从期初 t0 生存至某