周雪光事件史分析3

3．Cox模型 1 第三讲 COX模型：比例风险模型 本讲提要 1．Cox模型的内容 2．Cox模型的特点和应用 3．Cox模型的假设和检验 1. 什么是Cox模型 A。在事件史分析中的时间依赖问题 原因：时间是一个重要维度 例子： • 职务升迁：对工作期限的依赖 • 工作变动：对年龄（工龄）的依赖 三种解决办法: • 离散时间的事件史模型： o （1）时间作为解释变量，（2）使用随时间变化的变量，（3）不对时间 模型化 o 假设：时间依赖问题通过模型中的控制变量可以排除； • Cox模型：假设有时间依赖，但是把时间依赖看作是一个可以通过统计处理而绕 过的问题 • 参数方程：对时间依赖的方式加以模型化 B．COX模型 几个假设  假设我们讨论的事件史中，时间是连续的；  假设样本中每个人的风险是成比例的， 即每个人的风险都是其他人风险的一个固定比 例。所以，COX模型又称比例风险模型。 比例风险模型可以表述如下： h(t, x) = h0(t) eβ’x  ln h(t, x) = α(t) +β’x 3．Cox模型 2 在这里，α(t)=logh0(t)。 上面公式说明，一个样本在时间t的风险是两个部分的乘积： 1．h0(t) 是一个基准风险函数，它的形式没有被具体规定（不能是负数），可以用任 何形式出现。我们可以把它看作是当所有变量为0时，一个样本所面临的风险函 数。 2．eβ’x。一组k自变量的线形方程。 用广义线形模型（Generalized Linear Models）的语言来说，时间与解释变量间没有 interaction。.  If α(t) = α  exponential model  If α(t) = α t  Gompertz model  If α(t) = α log t  Weibull model 注意：Cox模型的一个重要假设是，h0(t)可以有任何形式，但是它对于所有的样本都是同 样的，因此，我们可以通过样本之间的相比来消除α(t)。这个假设的另外一个表述是，样 本之间的风险是成比例的（proportional）。 例如， 在这里，λ0(t)在公式中被取消，因此风险的比例在不同时间都是同样的，不随时间变 化。如果我们观察不同群体的log h(t)，那么它们的log h(t)应该是平行的。而h0(t) 可以 被解释为所有样本共有的共同基准风险。 logh(t) 3．Cox模型 3 2．Cox模型的特点和应用 （1） Cox的重要贡献 1．允许风险率有时间依赖，但对风险率的分布形式没有做任何假设。在h0(t)未知的 条件下，可以对上式中的β进行估算。 被称为部分似然值法。 2．统计计算上简易； 注意：部分似然值法的估算使用了有关持续期时间（duration times）的顺序，但对 事件发生的准确时间不予关注。事件而不是个人是似然值估算的关注点。 h0(t) 没有模型化，因此，Cox模型又被称为半参数模型。 估算方法被称为部分似然值法。对于“部分似然值法”，我们可以写为观察到的所有事件 的似然值的乘积。因此，如果J是事件的数目，那么我们可以有如下公式， 在这里，Lj是第j个事件的似然值。 （2） 参数估算方法（Estimation Method ） 在一个有n人的随机样本中，有k观察到的不同的生命时间（没有同时发生的事件）， n-k删截的生命时间。我们把这些k生命时间按时间长短从短到长顺序排列起来，即ti < tj if i < j. 设 Ri =R(ti) 为在时间ti的风险集，即在ti.之前仍然生存而且没有被删截的样 本。 当任何生命时间(ti) 出现时（即事件发生时），总是有一个相应的Ri. 个人 i [i ∈ Ri] 经 历该事件，而不是其他在风险集中的人，的概率是: 在这里，(h0(t))被取消了。 我们可以使用似然值法对下式求 β的估算值: 3．Cox模型 4 上式是所谓的“部分似然值法”，只对exp(xi'β)有关。而通常的似然值法需要知道 h0(t). 但在实际估算中，我们可以把上式看作为似然值函数，其β的估测具有最大似然 值法的各种特点： consistency and asymptotic normality。只是有着极少的效率损 失。 一个例子：STATA p22-23 id time x 1 1 3 2 5 2 3 9 4 4 20 9 5 22 10  h(subject 4 at time 20) = h4(20)=h0(20)exp(β0+9βx) 注意：在我们考虑 hj(20) 时，只有hj(20) 和hj(22)在风险集中。 我们可以用“部分似然值法”对上式求βx。 （3） 同时发生的事件（tied events） 3．Cox模型 5 如果时间是连续的，那么应该没一个以上的事件在同一时间发生。但是，在实际资料中， 时间并不是真正连续的。因此，当我们用一个时间整数来统计时间时，就会发生不同事件 在同一时间发生的情形。例如，可能几个人在一年里经历了“提升”事件。 Cox模型在部分似然值法的统计分析中，假设生命时间（ti）是各自不同的，在统计处理 时首先将时间T按顺序排列。我们可能碰到一个问题：如果资料中的样本经历的事件同时 发生，也就是说，在某一时间的T多次出现（tied events），怎么办？ 统计学家提出了几种近似方法来解决这个困难：  Breslow approximation  Marginal calculation （exactm）  Partial calculation （exactp）  Efron approximation （efron） （4） COX模型没有截距 Cox模型的统计结果中没有“截距”，这是因为截距是 α（基准风险函数）的一个部 分，在“部分似然值法”中被取消了。 （5）COX模型的问题 (1) Cox模型没有充分利用信息。只考虑持续时间的次序，但没有使用事件发生准确 时间的信息。 (2) 假设时间是连续的。对“同时发生的事件”处理困难。（但这些因素现在都不是 重要问题了。） (3) 没有直接计算 h0(t). （但是可以用其他方法估算出来。） 3．Cox模型的假设和检验 关于“比例风险”的假设的讨论 1．这个假设是什么意思？ 2．如何检验这个假设？ 3．Cox模型 6 我们从上式可以看出，比例风险意味着不随时间变化的相对风险，或者风险比率（hazard ratios）。对于一个定类变量 xj, exp(βj)可以解释为在模型中的类别与参考类别之间的相对 风险。对于一个连续性变量xj, exp(βjxj) 可以看做是xj每一单位变化（Δxj）导致的相对风 险。 例如，如果 xj 是一个虚拟变量（dummy variable）, xj = 1. 这表明，风险率的群体差异 是成比例的，即log(风险率) 的群体差异是平行的。时间与解释变量之间没有交互作用。 如果“比例风险”的假设不能成立，这相当于时间与某一个或几个变量相互作用。比例风 险模型假设每一个变量的作用在所有时间里都是相同的。如果一个变量的作用随时间而变 化，那么比例风险的假设就不能成立。 检验“比例风险”的假设 A．利用群体间log(风险率)平行的特性，画图检验，或使用log-rank检验等方 法。参见第一讲的“描述性统计方法”。 B．“再估测”检验：  思路：假设-使用的变量是适当的，检验其作用方式（functional form）  具体操作 Null hypotheses: β1 = 1 and β2 = 0. STATA output Cox regression -- Breslow method for ties No. of subjects = 494 Number of obs = 9932 No. of failures = 84 Time at risk = 9932 LR chi2(2) = 56.35 Log likelihood = -457.29343 Prob > chi2 = 0.0000 ----------------------------------------------------------------------------- _t | Coef. Std. Err. z P>|z| _hat | 1.787253 .5378555 3.32 0.001 _hatsq | -.1899061 .1270612 -1.49 0.135 ----------------------------------------------------------------------------- 3．Cox模型 7 C．将变量与时间相乘。在比例风险假设下（风险不随时间变化），这些相乘项应该没有 统计显著性；  例如：female*duration，educyr*duration D．Schoenfeld residuals 检验 基本思路：对COX模型估算后的 the residuals建立一个与时间有关的平滑函数，检验是否 有显著的统计关系。 (STATA，p. 160) 画图看：在“比例风险”的null假设下，我们期待这一函数的斜率为零（经过光滑处理 后）。 Null假设的检验： the log hazard ratio function 不随时间变化。 3．Cox模型 8 STATA output Cox regression -- Breslow method for ties No. of subjects = 494 Number of obs = 9932 No. of failures = 84 Time at risk = 9932 LR chi2(4) = 37.11 Log likelihood = -466.91272 Prob > chi2 = 0.0000 _t Haz. Ratio Std. Err. z P>z [95% Conf. Interval] rh female .1830379 .0849029 -3.66 0.000 .0737407 .4543332 educyr 1.11565 .0854919 1.43 0.153 .9600646 1.296448 t female 1.070113 .0317435 2.28 0.022 1.009671 1.134173 educyr 1.004917 .0049615 0.99 0.320 .9952394 1.014688 Note: Second equation contains variables that continuously vary with respect to time; variables are interacted with current values of _t. . test [t]; ( 1) [t]female = 0 ( 2) [t]educyr = 0 chi2( 2) = 6.28 Prob > chi2 = 0.0433 stphtest, detail; Test of proportional hazards assumption Time: Time ---------------------------------------------------------------- rho chi2 df Prob>chi2 ------------+--------------------------------------------------- female 0.26208 5.66 1 0.0173 educyr 0.13303 1.13 1 0.2876 ------------+--------------------------------------------------- global test 6.80 2 0.0333 ---------------------------------------------------------------- 3．Cox模型 9 3．Cox模型 10 3．Cox模型 11 资料设置和处理 统计分析 Cox regression -- Breslow method for ties No. of subjects = 494 Number of obs = 9932 No. of failures = 84 Time at risk = 9932 LR chi2(13) = 53.70 Log likelihood = -458.61947 Prob > chi2 = 0.0000 ------------------------------------------------------------------------------ _t | Haz. Ratio Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- female | .4475336 .1050049 -3.43 0.001 .2825582 .7088319 educyr | 1.156856 .0511697 3.29 0.001 1.060789 1.261622 govt | 5.242594 2.390388 3.63 0.000 2.145045 12.81315 public | 1.751634 .521873 1.88 0.060 .9768787 3.140842 firmcl | .8552261 .4045196 -0.33 0.741 .338427 2.16121 firmpr | .6173002 .3686971 -0.81 0.419 .191468 1.9902 firmfo | 1.547735 1.597446 0.42 0.672 .204718 11.70139 firmoth | 2.379482 1.731578 1.19 0.234 .571548 9.906312 fedmidl | .9897395 .2881935 -0.04 0.972 .5593286 1.751357 fedhigh | .9154186 .3243726 -0.25 0.803 .4570903 1.833317 fedcolg | .3554843 .1921114 -1.91 0.056 .1232573 1.025246 fedothr | .7052626 .3097568 -0.80 0.427 .2981924 1.668035 fparty | 1.711473 .4167235 2.21 0.027 1.06197 2.758212 ------------------------------------------------------------------------------ 实证假设检验 STATA output test govt=public; ( 1) govt - public = 0 chi2( 1) = 4.78 Prob > chi2 = 0.0289 . test fparty=2; ( 1) fparty = 2 chi2( 1) = 36.08 Prob > chi2 = 0.0000 3．Cox模型 12 不同的处理“tied events”的方法 非“比例风险”的Cox模型（nonproportional hazards） 如果发现“比例风险”的假设不能成立，有以下两种处理方法。 （1）分层（stratified）分析 在“分层”的情形下，我们将每个人的基准风险率相同的假设放松，让每个人的基准风险 率在自己的群体内相同，但是可以在不同群体间不同。 允许h(0)在男女之间不同，但是其他变量的测算不分男女。 stcox x1 x2 x3, strata(female) STATA output Stratified Cox regr. -- Breslow method for ties No. of subjects = 494 Number of obs = 9932 No. of failures = 84 Time at risk = 9932 LR chi2(12) = 40.48 Log likelihood = -405.06742 Prob > chi2 = 0.0001 ------------------------------------------------------------------------------ _t | Haz. Ratio Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- educyr | 1.155365 .0513153 3.25 0.001 1.059042 1.260449 govt | 5.557202 2.553305 3.73 0.000 2.258209 13.67566 public | 1.830654 .5477345 2.02 0.043 1.018417 3.290691 firmcl | .8408571 .3985934 -0.37 0.715 .3320666 2.129213 firmpr | .6661351 .3985283 -0.68 0.497 .2062122 2.151841 firmfo | 1.890063 1.964058 0.61 0.540 .2465761 14.48778 firmoth | 2.276913 1.660455 1.13 0.259 .5452571 9.508048 fedmidl | .9752392 .2842765 -0.09 0.931 .5507962 1.726758 fedhigh | .8821322 .3130408 -0.35 0.724 .4400167 1.768472 fedcolg | .3464516 .1875487 -1.96 0.050 .1199089 1.000999 fedothr | .6967778 .3071841 -0.82 0.413 .2936502 1.653325 fparty | 1.68884 .413863 2.14 0.032 1.044708 2.730121 ------------------------------------------------------------------------------ Stratified by female （2）与时间的相交作用（interaction） 3．Cox模型 13 统计结果的解释 定类变量（性别、工作单位）：注意“参考群体” 连续变量 log relative hazard (LRH): exp(β) 其他有关模型测算的统计分析数据 -log-L χ2 COX模型在离散性时间模型中的推广 Logistic模型在严格意义上不是比例风险模型。但当时间区间很小，风险很小的条件下， logit模型近似“比例风险”模型。 如果假设比例风险模型描述了观察到的事件发生的过程，那么我们可以从统计理论上推 论，cloglog模型是Cox模型应用在离散性时间的相应模型。 3．Cox模型 14 作业3：COX模型的统计计算和解释 继续使用你在作业2中的资料和变量。应用COX模型对同一资料和因果关系进行统 计分析。写一个简要的报告： 1．讨论COX模型对你的研究的适用性； 2．讨论统计模型解释能力、实证假设检验、变量参数解释； 3．检验COX模型有关比例风险的假设； 4．比较离散时间的logistic模型和Cox模型分析的结果。 3．Cox模型 15 /*****************************************************************************/ /* */ /* EHA_COX.DO */ /* */ /* This program illustrates the estimation of the Cox model, and various */ /* model building and diagonistic techniques. */ /* 01/05/05 xz; */ /* */ /*****************************************************************************/ #delimit ; clear; set mem 50m; set more off; set matsize 800; set more off; * First, link to the dataset, and define working directory; * use c:\data\sample600_eha, clear; global log c:\data; global data c:\data; cap log close; log using $log\eha_cox.log, replace; use $data\eha600_new; /*****************************************************************************/ /* */ /* The program below uses STATA's 'stset' procedure to define the event */ /* history dataset. Once the event history data are properly defined in */ /* STATA, you canconduct various statistical analyses using the dataset. */ /* */ /*****************************************************************************/ /* An illustration using party membership as dependent var */ * Define a variable 'riskyear' for the year when one is at risk: age=18; * select appropriate riskset: age >= 18; keep if t0 - birth >= 18; * We now estimate the Cox model; * Note: unlike other statistical models estimated using STATA, once you stset your event history data, you do not need to specify the 'response variable' in the model specification below. This also means that you need to 'stset'(that is, to specify new 3．Cox模型 16 response variable and the data structure your data whenever you analyze a different response variable; * The final model; stcox female educyr govt public firmcl firmpr firmfo firmoth fedmidl fedhigh fedcolg fedothr fparty; test govt=public; test fparty=2; * Using different tie-handling approximation: Efron (Default: breslow. Other options are: exactm exactp.); stcox female educyr govt public firmcl firmpr firmfo firmoth fedmidl fedhigh fedcolg fedothr fparty, efron; * Obtaining stratified analysis; stcox educyr govt public firmcl firmpr firmfo firmoth fedmidl fedhigh fedcolg fedothr fparty, strata(female); * Obtaining estimates of baseline functions--This program does not work!!!; *stcox period educyr, strata(female) basech(H0); *graph H0 _t, c(J) sort; *gen H0 = H if female == 0; *gen H1 = H if female == 1; *graph H0 H1 _t, c(JJ) sort ylabel xlabel; * Diagnosis of Cox model; * First, estimate the model; stcox female educyr govt public firmcl firmpr firmfo firmoth fedmidl fedhigh fedcolg fedothr fparty; * Tests based on re-estimation; linktest; * Test based on interaction with analysis time; * According to STATA, 'tvc(x) specifies that I want x to be interacted with some function of analysis time,' 'texp(_t) specifies the function to be used'; stcox female educyr, tvc(female educyr) texp(_t); test [t]; *test based on Schoenfeld residuals (I use a simplified model); 3．Cox模型 17 stcox female educyr, schoenfeld (sch*) scaledsch(sca*); stphtest, detail; * Graphical methods; stphtest, plot(educyr) title("PH test based on Schoenfeld residtual"); gr save stphtest_educyr, replace ; stphtest, plot(female) title("PH test based on Schoenfeld residtual"); gr save stphtest_female, replace; stphplot, by(female) title("Graph by female to check for parallel lines"); gr save stph_sex, replace; stphplot, strata(female) adjust(educyr) title("Graph by female for parallel lines, controling for educyear"); gr save stph_sex_ed, replace; 3．Cox模型 18 --------------------------------------------------------------------------------------------- --- log: c:\data\eha_cox.log log type: text opened on: 18 Jan 2005, 07:31:23 . use $data\eha600_new; . /*****************************************************************************/ > /* */ > /* The program below uses STATA's 'stset' procedure to define the event */ > /* history dataset. Once the event history data are properly defined in */ > /* STATA, you canconduct various statistical analyses using the dataset. */ > /* */ > /*****************************************************************************/ > > /* An illustration using party membership as dependent var */ > > * Define a variable 'riskyear' for the year when one is at risk: age=18; . * select appropriate riskset: age >= 18; . keep if t0 - birth >= 18; (600 observations deleted) . * We now estimate the Cox model; . * Note: unlike other statistical models estimated using STATA, > once you stset your event history data, you do not need to > specify the 'response variable' in the model specification below. > This also means that you need to 'stset'(that is, to specify new > response variable and the data structure your data whenever you > analyze a different response variable; . * The final model; . stcox female educyr > govt public firmcl firmpr firmfo firmoth > fedmidl fedhigh fedcolg fedothr fparty; failure _d: party analysis time _t: (year-origin) origin: time riskyear id: id Iteration 0: log likelihood = -485.47001 Iteration 1: log likelihood = -464.10861 Iteration 2: log likelihood = -458.69364 Iteration 3: log likelihood = -458.61955 Iteration 4: log likelihood = -458.61947 Refining estimates: Iteration 0: log likelihood = -458.61947 Cox regression -- Breslow method for ties 3．Cox模型 19 No. of subjects = 494 Number of obs = 9932 No. of failures = 84 Time at risk = 9932 LR chi2(13) = 53.70 Log likelihood = -458.61947 Prob > chi2 = 0.0000 ------------------------------------------------------------------------------ _t | Haz. Ratio Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- female | .4475336 .1050049 -3.43 0.001 .2825582 .7088319 educyr | 1.156856 .0511697 3.29 0.001 1.060789 1.261622 govt | 5.242594 2.390388 3.63 0.000 2.145045 12.81315 public | 1.751634 .521873 1.88 0.060 .9768787 3.140842 firmcl | .8552261 .4045196 -0.33 0.741 .338427 2.16121 firmpr | .6173002 .3686971 -0.81 0.419 .191468 1.9902 firmfo | 1.547735 1.597446 0.42 0.672 .204718 11.70139 firmoth | 2.379482 1.731578 1.19 0.234 .571548 9.906312 fedmidl | .9897395 .2881935 -0.04 0.972 .5593286 1.751357 fedhigh | .9154186 .3243726 -0.25 0.803 .4570903 1.833317 fedcolg | .3554843 .1921114 -1.91 0.056 .1232573 1.025246 fedothr | .7052626 .3097568 -0.80 0.427 .2981924 1.668035 fparty | 1.711473 .4167235 2.21 0.027 1.06197 2.758212 ------------------------------------------------------------------------------ . test govt=public; ( 1) govt - public = 0 chi2( 1) = 4.78 Prob > chi2 = 0.0289 . test fparty=2; ( 1) fparty = 2 chi2( 1) = 36.08 Prob > chi2 = 0.0000 . * Using different tie-handling approximation: Efron > (Default: breslow. Other options are: exactm exactp.); . stcox female educyr > govt public firmcl firmpr firmfo firmoth > fedmidl fedhigh fedcolg fedothr fparty, efron; failure _d: party analysis time _t: (year-origin) origin: time riskyear id: id 3．Cox模型 20 Iteration 0: log likelihood = -485.06059 Iteration 1: log likelihood = -463.57643 Iteration 2: log likelihood = -457.98454 Iteration 3: log likelihood = -457.89975 Iteration 4: log likelihood = -457.89965 Refining estimates: Iteration 0: log likelihood = -457.89965 Cox regression -- Efron method for ties No. of subjects = 494 Number of obs = 9932 No. of failure