Stata: 生存分析一文读懂

发布时间：2020-03-03 阅读 9688

Stata 连享会主页 || 视频 || 推文

温馨提示： 定期清理浏览器缓存，可以获得最佳浏览体验。

✌ 课程详情： https://gitee.com/arlionn/Course | lianxh.cn

✨ 课程主页： https://gitee.com/arlionn/Course

作者：李琼琼 (山东大学)
Email: lqqflora@163.com

1. 生存分析简介

生存分析法 (survival analysis) 源于生物统计学和人口学，主要对人或动物的寿命或生存时间进行分析，现在也广泛应用于社会科学和商业等领域。Cameron and Trival (2005) 从计量的角度认为生存分析主要研究 从一种状态转变到另一种状态所经历的时间。其应用非常广泛，比如：

失业持续的时间
初创企业存活时间
从开始筹备到顺利 IPO 所需的时间

目前，国内已有较多的学者使用生存分析来研究中国出现的社会经济问题,研究主题有：

农村劳动力向城市流动的影响因素分析（张世伟, 赵亮, 2009）
企业生存分析 (吴斌, 2006; 逯宇铎等，2013; 许家云，毛其淋，2016; 刘海洋等, 2017)
扭曲工资与企业成长 (吴先明, 2017)
企业新产品创新 (毛其淋,许加云, 2014; 毛其淋, 许加云, 2015)
企业出口 (毛其淋, 盛斌, 2013)
中国进出口贸易持续期 (邵军, 2011; 谭晶荣, 童晓乐, 2014; 何树全, 张秀霞, 2011; 赵瑞丽等, 2016)
上市公司 "ST特别处理" (吕长, 江赵岩, 2004)
企业财务困境预警 (宋雪枫,杨朝军,2006; 陆志明等，2007; 王晓鹏等, 2007; 陆志明等, 2007; 过新伟, 胡晓, 2012; 李宏伟,2019)
网络贷款违约风险 (李思瑶等, 2016)
女性生育间隔 (靳永爱等, 2019)
……

1.1 基本介绍

生存分析 (survival analysis) 又被称为久期分析 (duration analysis) 、转换分析 (transition analysis)、风险分析 (hazrd analysis)、信度分析(reliability analysis)、失效时间分析(failure-time analysis)、历史事件分析 (event history analysis) 等，主要研究 关注事件 发生所需要的时间，具体可归纳为三类 (Paul D. Allision,2018)：

定性变量状态变化所花费的时间, 比如离婚、晋升、死亡
定量变量出现急剧变化所经历的时间, 比如人口总和出生率锐减 (见 Figure 1)
定量变量超过某个阈值所需要的时间, 比如将体重减到健康标准 (见 Figure 2)

Figure 1: Quantitative variables change sharply

Figure2: Quantitative variables cross a threshold

1.2 生存分析解决的问题

例子

有 400 名犯人刑满被释放，从监狱出来后，研究人员开始追踪调查他们的情况，追踪时间为 2 年，以期观察他们是否会再次被捕。在这里，研究人员关注事件是 再次被逮捕，解释变量为救助金、教育、工作状态等。

研究方法 A：用二元被解释变量 (逮捕 = 1, 没有逮捕 = 0) 对解释变量进行线性回归；
- 局限：应该使用逻辑回归 (Logistic Model); 没有充分利用信息; 不能有效处理工作状态这一变量(此变量随时间而变动)
研究方法 B：把被解释变量定义为犯人从释放到第一次被捕的时间，再进行线性回归。
- 局限：出现有些个体没有被逮捕或者在追踪期以后被逮捕却被记录成 2 年的情况 (censoring), 并且依旧无法概括随时间变动的解释变量。

评价: 对于 "关注事件发生所需要的时间" 这一类的研究，用传统模型做实证，会面临截堵 (censoring), 以及解释变量发生了时变等主要问题, 而使用生存分析法则可以解决这些问题。

截堵 (censoring)

截堵分为 3 种类型：即右截堵 (right censoring)、左截堵 (left censoring) 和区间截堵 (interval censoring), 下面以生存分析的例子来进行区分。

设定一个恒定的值 $c$ (被解释变量，持续的时间)

右截堵: 将 $c$ 作为被解释变量的上限，但是存在一些观测值 $T$ 大于 $c$ 。例如: 对一组年龄为 30 岁的女性进行采访, 关注的事件是初婚时间。采访中提出的问题是「结婚时的年龄 ( $T$ )」。本例中 $c = 30$ 。然而，若样本中包含部分未婚女性, 她们的结婚年龄一定会超过 30 岁, 但由于抽样方法的局限，其具体结婚年龄无法在访谈中得知。
左截堵: 将 $c$ 作为被解释变量的下限，但是存在一些观测值 $T$ 小于 $c$ 。例如: 对一组年龄在 20 岁及以上的女性进行采访, 关注的事件同样是初婚年龄。最开始提的问题为是否结婚, 这里的 $c$ 为 20。但是有一些 20 岁的女性已经结婚，则她们的结婚年龄一定小于 20, 具体是多少无法从提问中知晓。
区间截堵: 将 $c_{1}$ 和 $c_{2}$ 分别作为被解释变量的下限和上限, 存在观测值 $T$ 小于 $c_{1}$ , 也存在观测值 $T$ 大于 $c_{2}$ 。例如: 对一组年龄在 20 岁到 30 岁之间中国女性进行采访，问她们第一次生育的年龄, 这里 $c_{1}$ 等于 20 (中国女性合法结婚年龄, 婚后生孩子才是合法的)， $c_{2}$ 为 30。 30 岁的女性还没有生育，则其生育年龄一定超过 30, 对于 20 岁以前生孩子的女性，考虑到隐私及法律问题，并没有回答具体生育年龄, 只是说了已经有孩子, 则她们的生育年龄一定小于 20。

其中，右截堵是最常见的类型, 也是生存分析主要研究对象。

2. 生存分析模型

2.1 基本概念

假设: 在生存分析中, 持续时间 $T$ 是随机的。

累积分布函数

F (t) = \Pr (T \leq t)

含义: 在时间 $t$ 或 $t$ 以内，失败 (死亡) 的概率。

由此可以定义关注事件发生的概率对应的密度函数: $f (t) = d F (t) / d t$

生存函数

S (t) = 1 - F (t)

含义: 生存超过某个时间 $t$ 的概率，即关注事件在时间 $t$ 或 $t$ 以内都未发生的概率

风险函数

风险函数 $λ (t)$ 关注事件在 $t$ 时刻发生的概率, 也被称为在 $t$ 时刻 (瞬间) 发生的风险。推导过程如下：

累积风险函数

\begin{aligned} Λ (t) & = \int_{0}^{t} λ (s) d s \\ = - \ln S (t) \end{aligned}

含义: 关注事件到 $t$ 时刻为止发生的概率，相比较风险函数更容易被精确估计。

2.2 风险函数的类型

指数模型 (Exponential model)

当 $T$ 随机变量服从指数分布, 参数为 $λ$

λ (t) = λ

解释: 指数分布的风险函数为常数, 说明瞬间关注事件发生的概率并不依赖已持续的时间。

冈珀茨模型 (Gompertz model)

当 $T$ 随机变量服从冈珀茨分布, 参数为 $λ$ 和 $α$

λ (t) = λ_{0} e^{α t} where λ_{0} = e^{μ}

解释：冈珀茨分布的风险函数呈指数变化率 $e^{α t}$ , 并为单调函数:

当 $α < 0$ 时，单调递减;
当 $α = 0$ 时，变为指数模型的风险函数;
当 $α > 0$ 时，单调递增。

威布尔模型 (Weibull model)

当 $T$ 随机变量服从威布尔分布, 参数为 $λ$ 和 $α$

λ (t) = λ_{0} t^{α} where α > - 1

解释：与冈珀茨模型不同, 威布尔分布的风险函数呈幂函数变化率 $t^{α}$ , 也为单调函数:

当 $α < 0$ 时，单调递减;
当 $α = 0$ 时，变为指数模型的风险函数;
当 $α > 0$ 时，单调递增。

小结

以上 3 种模型统一被称为 比例风险模型 (proportional hazards models, PH), 是随机变量 $T$ 符合特定假设的参数模型, 均使用最大似然法估计 (MLE)。其对数形式如下:

指数模型： $\log λ (t) = β_{0} + β_{1} x_{1} + β_{2} x_{2}$
冈珀茨模型： $\log λ (t) = β_{0} + β_{1} x_{1} + β_{2} x_{2} + α t$
威布尔模型： $\log λ (t) = β_{0} + β_{1} x_{1} + β_{2} x_{2} + α \log t$

注:

$β$ 为半弹性, 即某个解释变量增加 1 个单位，将导致风险函数平均增加百分之 $β$ 。
$e^{β}$ 为风险比率 (Hazard Ratio), 即某个解释变量变量增加 1 单位, 将导致新风险比率变为原来的 $e^{β}$ 倍。
另外, 不含解释变量 $X$ 的风险函数通常被称为 基准风险 $λ_{0} (t)$ 。

此外，如果使用比例风险模型，需要满足 比例风险假定, 即风险函数的形式为:

λ (t; x) = λ_{0} (t) h (x) = λ_{0} (t) e^{x^{'} β}

检验的方式有: 对数-对数图 (log-log plot), 观测-预测图 (observed versus expected plot), 舍恩菲尔德残差检验 (Schoenfeld residuals tests)

2.3 加速失效时间模型

加速失效时间模型 (accelerated failure-time model, AFT) 主要研究解释变量 $x$ 对平均寿命 (从关注事件未发生到发生平均经历的时间 $T$ ) 的影响, 而上节提到的比例风险模型分析的角度则为解释变量 $x$ 对风险函数 $λ (t; x)$ 的作用。

加速失效时间模型的设定为:

\ln t = x^{'} β + u

其中, $u$ 分布不同，会形成不同的 AFT 模型。 $b e t a$ 的含义和上一节也不相同，在这里 $b e t a$ 表示某个解释变量 $x$ 增加一个单位, 能使平均寿命增加百分之 $b e t a$ 。

为何被称为 "加速失效时间模型" ?

持续时间 (生存时间):

t = \exp (x^{'} β) v, where v = e^{u}

那么，

v = t \cdot \exp (- x^{'} β)

风险函数:

λ (t; x) = λ_{0} (v) \exp (x^{'} β)

把 $v$ 代入后

λ (t; x) = λ_{0} (t \exp (- x^{'} β)) \exp (x^{'} β)

风险函数满足上面的形式，则为加速失效时间模型。若 $\exp (- x^{'} β) > 1$ , 称为加速， $\exp (- x^{'} β) < 1$ 则为减速，不过无论 "加速" 还是 "减速"，都是加速失效时间模型。

哪些具体的模型可以视为加速失效时间模型 ?

Cameron and Trival (2005) 总结有 5 种加速失效时间模型, 其中指数模型和威布尔模型既是比例风险模型 (PH) 又是加速失效时间模型 (AFT).

![Table1：参数模型](x https://fig-lianxh.oss-cn-shenzhen.aliyuncs.com/Survival-Analysis-AFT.png)

Table 17.5. Standard Parametric Models and Their Hazard and Survivor Functions

Parametric Model	Hazard Function	Survivor Function	Type
Exponential	$γ$	$\exp (- γ t)$	PH, AFT
Weibull	$γ α t^{α - 1}$	$\exp (- γ t^{α})$	PH, AFT
Generalized Weibull	$γ α t^{α - 1} S (t)^{- μ}$	${[1 - μ γ t^{α}]}^{1 / μ}$	PH
Gompertz	$γ \exp (α t)$	$\exp (- (γ / α) (e^{α t} - 1))$	PH
Log-normal	$\frac{\exp (- (\ln t - μ)^{2} / 2 σ^{2})}{t σ \sqrt{2 π} [1 - Φ (\ln t - μ) / σ)]}$	$1 - Φ ((\ln t - μ) / σ)$	AFT
Log-logistic	$α γ^{α} t^{α - 1} / [(1 + (γ t)^{α})]$	$1 / [1 + (γ t)^{α}]$	AFT
Gamma	$\frac{γ (γ t)^{α - 1} \exp [- (γ t)]}{Γ (α) (1 - I (α, γ t)]}$	$1 - I (α, γ t)$	AFT

Source: Cameron and Trival (2005), Table 17.5

2.4 Cox PH 模型

比例风险模型和加速失效时间模型均为参数模型, 需要对风险函数的具体形式作出假设，再用最大似然法估计。但是截堵数据可能会导致风险函数设定错误, 此时会出现不一致的 MLE 估计。Cox (1972, 1975) 以 PH 模型为基础提出估计 $β$ 的半参数模型, 被定义为 Cox PH 模型或者 Cox 模型。Cox PH 模型不需要假设基准风险 $λ_{0} (t)$ 的具体形式，因为个体 $i$ 和个体 $j$ 的风险函数之比可以将 $λ_{0} (t)$ 约去, 只剩下与解释变量 $x$ 相关项。

\frac{λ (t; x_{i})}{λ (t; x_{j})} = \frac{λ_{0} (t) e^{x_{i}^{i} β}}{λ_{0} (t) e^{x_{j}^{'} β}} = e^{{(x_{i} - x_{j})}^{'} β}

由于构成风险函数 $λ_{0} (t) e^{x^{'} β}$ 前部分 $λ_{0} (t)$ 不需要估计参数 (非参数部分) 而后半部分 $e^{x^{'} β}$ 需要估计参数，故 Cox PH 模型是一种半参数回归，同时要和 PH 模型一样, 风险函数满足 $λ (t; x) = λ_{0} (t) e^{x^{'} β}$ 的设定形式。

此外，在过去的推文曾经介绍过使用 Tobit 模型来解决数据截堵的问题, Cox PH 模型和 Tobit 模型的方法存在很大差异，Tobit 模型无法很好地估计风险率，而且总体上学者们对 Cox PH 模型的认可度更高 (Cameron and Trival, 2005)。

3. 生存分析的基本步骤及 Stata 命令

本章节主要介绍做生存分析的基本步骤及实现的 Stata 命令, 以便于更好地分析生存数据。在介绍完操作程序后, 再通过介绍中国工业经济的一篇文章, 来让大家对生存分析的应用有更好的理解。基本步骤有：

第一步 设定生存分析数据
第二步 画生存函数、累积风险函数和风险函数
第三步 进行参数回归及非参数回归
第四步 比例风险假定的检验

使用的数据是来自 Stata-press 官网的 hip.dta 数据（右击可以另存到本地）。该数据集主要调查一种新的可穿戴充气防护装置是否可以保护老年人，即降低老年人发生跌倒骨折的风险。邀请了 48 位年龄为 60 岁以上的女性展开一项试验，任意挑选 28 名女性发放并让她们穿戴防护装置，作为实验组, 剩下的 20 名女性则不穿防护装置作为控制组。

研究中的关键变量定义如下：

关注事件 为骨折 (fracture)
持续时间 为开始被调查到发生骨折(fracture = 1) 的时间, 同时也包括发生过骨折以后恢复了再次发生骨折的时间 (被试会被调查几个阶段，不同阶段之间会有时间间隔)
解释变量 包括：年龄 (age)、血钙（calcium) 和穿戴防护装置 (protect)。

3.1 将数据设定为生存分析格式

首先，我们了解一下原始数据的概况：

. use http://www.stata-press.com/data/cggm3/hip, clear //调用 hip.dta 数据
. describe //对数据集进行描述
*-----------使用 describe 命令 table2 ------
Contains data from http://www.stata-press.com/data/cggm3/hip.dta
  obs:           106                          hip fracture study
 vars:             7                          30 Jan 2009 11:58
 size:         1,060                          
-------------------------------------------------------------------
              storage   display    value
variable name   type    format     label      variable label
-------------------------------------------------------------------
id              byte    %4.0g                 patient ID
time0           byte    %5.0g                 begin of span
time1           byte    %5.0g                 end of span
fracture        byte    %8.0g                 fracture event
protect         byte    %8.0g                 wears device
age             byte    %4.0g                 age at enrollment
calcium         float   %8.0g                 blood calcium level
-------------------------------------------------------------------

stset 命令

stset 命令将数据设定为适合生存分析的形式, 其语法为，

. stset tvar, failure(fail) id(idvar) origin(time torig)

各要素的含义如下：

tvar 为结束时间
failure(fail) 设定关注事件, 默认 fail==1 为关注事件发生
id(idvar) 设定样本中的用以标示每个观察值的变量
origin(time torig) 设定 torig 为开始时间, 默认为 0。

进行 stset 设定后，Stata 呈现结果如下：

. stset time1, origin(time time0) id(id) failure(fracture == 1)
*------------使用 stset 命令 table3 -------

                id:  id
     failure event:  fracture == 1
obs. time interval:  (time1[_n-1], time1]
 exit on or before:  failure
    t for analysis:  (time-origin)
            origin:  time time0

--------------------------------------------------------------------------
        106  total observations
          0  exclusions
--------------------------------------------------------------------------
        106  observations remaining, representing
         48  subjects
         31  failures in single-failure-per-subject data
        744  total analysis time at risk and under observation
                                                at risk from t =         0
                                     earliest observed entry t =         0
                                          last observed exit t =        39

据表 3 显示, 样本中共有 106 个观测值, 48 个被试, 在调查期间有 31 个被试发生过摔倒骨折。另外, 执行 stset 命令后会自动生成 4 个变量，分别为：_st, _d, _t, _t0，含义为:

_st = 1 表示该观测值的数据会被使用,
_d 表示在一个记录中是否发生过关注事件,
_t 代表一个记录的结束时间,
_t0 代表一个记录的开始时间。

stdescribe 命令

stdescribe 命令是在 stset 命令之后使用的, 用于描述生存分析数据的基本特征

. stdescribe
*------------使用 stdescribe 命令 table4 -------
         failure _d:  fracture == 1
   analysis time _t:  (time1-origin)
             origin:  time time0
                 id:  id

                                   |-------------- per subject --------------|
Category                   total        mean         min     median        max
------------------------------------------------------------------------------
no. of subjects               48   
no. of records               106    2.208333           1          2          3

(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

failures                      31    .6458333           0          1          1
------------------------------------------------------------------------------

据表 4 显示, 48 个被试一共被记录过 106 次, 平均记录超过 2 次, 一共骨折了 31 次, 平均每人发生 0.64 次骨折。第 3、4 行表示若记录开始平均为 0, 终止 (发生骨折,或者没有发生骨折但是退出记录) 时间平均为 15.5 个月, 即大概 48 个被试平均 15.5 个月会摔倒一次，最短为 1 个月，最长为 39 个月。

3.2 描述性分析 (非参数分析)

生存分析的第二步是画生存函数、累积风险函数和风险函数图，主要使用 sts graph 命令，可以绘制 Kaplan-Meier 生存估计量, Nelson-Aalen 累积风险估计量和风险函数图。

生存函数图

Kaplan-Meier 估计量大致等于样本存活 (关注事件未发生) 时间超过 $t$ 的观测值数目占样本总体观测值的数目比值, 此方法无需对数据的分布做假设、也不需要进行参数估计, 故被称为非参数分析。

*-----使用 sts graph 命令 figure5 -----
//数据已设置成生存分析数据
. use http://www.stata-press.com/data/cggm3/hip2, clear 
. sts graph  //生存分析图
. graph export figure5.png

从总体生存函数图无法得知防护装置对生存的功效，下面画分组的生存函数图。

*-----使用 sts graph 命令 figure6-----
* 以是否穿防护装置分组画生存分析图  
* plot2opts(lp("-"))表示第二个图用虚线表示
. sts graph, by(protect) plot2opts(lp("-")) 
. graph export figure6.png

从图 6 可以看出是否穿戴防护装置对生存函数有明显影响, 相比较控制组，实验组有更好的生存经历。

累积风险函数图

Nelson-Aalen 累积风险估计量是对局部风险率加总得出的，估计量为

\hat{Λ} (t) = \sum_{j | t_{j} \leq t} {\hat{λ}}_{j} = \sum_{j | t_{j} \leq r} \frac{d_{j}}{r_{j}}

*-----使用 sts graph 命令 figure7-----
sts graph, cumhaz ci //cumhaz 选项表示画累积风险函数, ci 显示 95% 的置信区间。
graph export figure7.png

风险函数图

风险率估计量可表示为:

λ_{j} = \frac{d_{j}}{r_{j}}

其中, $d_{j}$ 表示在时刻 $t_{j}$ 出现关注事件的个体的数量, $r_{j}$ 表示在时刻 $t_{j - 1}$ 还未出现关注事件 ("存活") 的个体的数量。

*-----使用 sts graph 命令 figure8-----
sts graph, cumhaz ci //cumhaz by(protect)选项表示画累积风险函数, ci 显示 95% 的置信区间。
graph export figure7.png

由图 8 可知，是否穿戴防护装置对风险函数影响较大, 控制组的发生骨折的风险一直远远高于实验组。

上述分析均没有考虑解释变量，只是对被解释变量进行初步的描述性分析, 下面进行更为准确的模型估计分析。

3.3 模型估计

比例风险模型

. streg protect age calcium, nohr nolog dist(weib)
*- dist(weib) 代表威布尔回归的风险函数, 
*- dist(e) 为指数回归，dist(gom) 为刚珀茨回归，
*- nohr 显示回归系数而不显示风险比率

*-----使用 streg 命令 table5-----
         failure _d:  fracture
   analysis time _t:  time1
                 id:  id

Weibull PH regression

No. of subjects =           48                  Number of obs    =         106
No. of failures =           31
Time at risk    =          714
                                                LR chi2(3)       =       35.22
Log likelihood  =   -41.687862                  Prob > chi2      =      0.0000

------------------------------------------------------------------------------
          _t |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
     protect |  -2.192387   .4078599    -5.38   0.000    -2.991778   -1.392996
         age |   .0801079   .0550936     1.45   0.146    -.0278735    .1880893
     calcium |  -.1744454   .2194006    -0.80   0.427    -.6044626    .2555718
       _cons |  -7.787834   5.734033    -1.36   0.174    -19.02633    3.450665
-------------+----------------------------------------------------------------
       /ln_p |   .5207567   .1369941     3.80   0.000     .2522531    .7892602
-------------+----------------------------------------------------------------
           p |   1.683301   .2306023                      1.286922    2.201767
         1/p |   .5940709   .0813842                      .4541807     .777048
------------------------------------------------------------------------------

解读: 据表 5 显示, 穿戴防护装置会导致老年人骨折的风险降低 2.19%, 在 1% 的显著性水平下显著, 年龄上升会导致骨折风险的增加, 血钙的增加会降低骨折, 但不显著。倒数第 3 行 ln_p = 0 的 p 值为 0, 显著拒绝指数回归的原假设, 认为应该使用威布尔回归。倒数第 2 行的 p 值大于 1, 说明风险函数随时间而递增。

*-----使用 stcurve 命令 figure 9-----
. stcurve, hazard //画风险函数图

加速时间失效模型 (AFT)

选择加速失效时间模型的对数正态回归

*-----使用 streg 命令 table6-----
.  streg protect age calcium,  nolog dist(logn)

         failure _d:  fracture
   analysis time _t:  time1
                 id:  id

Lognormal AFT regression

No. of subjects =           48                  Number of obs    =         106
No. of failures =           31
Time at risk    =          714
                                                LR chi2(3)       =       35.50
Log likelihood  =   -41.758706                  Prob > chi2      =      0.0000

------------------------------------------------------------------------------
          _t |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
     protect |    1.44531   .2495669     5.79   0.000     .9561681    1.934452
         age |  -.0639809    .041565    -1.54   0.124    -.1454468    .0174849
     calcium |   .0692953   .1668585     0.42   0.678    -.2577412    .3963319
       _cons |    5.74614   4.416337     1.30   0.193    -2.909722      14.402
-------------+----------------------------------------------------------------
    /lnsigma |   -.294334   .1271106    -2.32   0.021    -.5434661   -.0452018
-------------+----------------------------------------------------------------
       sigma |   .7450276   .0947009                      .5807319    .9558046
------------------------------------------------------------------------------

解读: 对数回归的系数符号与比例风险模型的系数符号相反, 在 2.3 小节曾做过解释, 比例风险模型研究风险率的变化而 AFT 研究的是平均持续时间。防护装置可以显著降低风险率, 增加平均持续时间。

*-----使用 stcurve 命令 figure 10-----
. stcurve, hazard //画风险函数图

Cox PH 模型

对参数回归的具体分布无法确定, 故下面使用半参数 Cox PH 回归。

. stcox protect age calcium,r nohr nolog
*-----使用 streg 命令 table 7-----
Cox regression -- Breslow method for ties

No. of subjects      =           48             Number of obs    =         106
No. of failures      =           31
Time at risk         =          714
                                                Wald chi2(3)     =       31.21
Log pseudolikelihood =   -82.062396             Prob > chi2      =      0.0000

                                    (Std. Err. adjusted for 48 clusters in id)
------------------------------------------------------------------------------
             |               Robust
          _t |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
     protect |  -2.242575   .4328035    -5.18   0.000    -3.090855   -1.394296
         age |   .0676642   .0438904     1.54   0.123    -.0183594    .1536878
     calcium |  -.2114179   .1600249    -1.32   0.186    -.5250609    .1022251
------------------------------------------------------------------------------

解读: 从表 7 可以看出, Cox PH 回归和前面参数回归的结果类似, 但不依赖于具体的分布假设, 故结果更稳健。

*-----使用 stcurve 命令 figure 11-----
. stcurve, hazard //画风险函数图

本文从介绍的角度讨论了以上模型，实际做实证时挑选适合的部分模型做基本模型即可, 挑选的规则有对数似然值和 AIC 等。

3.4 模型的检验

比例风险模型和 Cox PH 模型的重要假设是 $λ (t; x) = λ_{0} (t) e^{x^{'} β}$ , 如果这个假设不成立则上述模型不成立。因此，必须要对比例风险模型设定进行检验。下面介绍三种检验方式。

对数 - 对数图

若对数 - 对数图中的曲线相互平行, 则比例风险设定是成立的。

*-----使用 stphplot 命令 figure 12-----
. stphplot, by(protect)

两条曲线大致是平行的, 满足比例风险假定。

观测 - 预测图

若观测值和预测值非常接近, 则比例风险设定是成立的。

*-----使用 stcoxkm 命令 figure 13-----
. stcoxkm, by(protect)

生存函数的预测值和观测值总体上相距比较近, 基本满足比例风险假定。

舍恩尔德残差图

对于检验比例风险假定最有用的是舍恩尔德残差，如果比例风险假设成立, 则舍恩尔德残差不应该随时间呈现规律性变化。对于每个变量 $x$ 都可以将其舍恩尔德残差对时间回归, 并考察其斜率是否为 0, 也可以将其舍恩尔德残差和时间画图, 观察其斜率是否为 0.

quitely stcox protect age calcium,r nohr nolog
estat phtest, detail

*-----使用 estat phtest命令 table8 -----    
    Test of proportional-hazards assumption

    Time:  Time
    ----------------------------------------------------------------
                |       rho            chi2       df       Prob>chi2
    ------------+---------------------------------------------------
    protect     |      0.00064         0.00        1         0.9974
    age         |     -0.10011         0.17        1         0.6828
    calcium     |      0.00017         0.00        1         0.9997
    ------------+---------------------------------------------------
    global test |                      0.26        3         0.9674
    ----------------------------------------------------------------

三个解释变量都没有拒绝原舍恩尔德残差对时间回归的斜率为 0 的原假设, 故支持比例风险假定。

*-----使用 estat phtest命令 figure14 -----   
. estat phtest, plot(protect)

从残差与时间的拟合图来看其斜率大致为 0，满足比例风险假定。

3.5 国内研究实例

吴先明等 (2017) 考察了工资扭曲 (劳动者实际所得收入小于劳动价值应得收入) 与企业成长之间的关系, 并引入种群密度研究三者之间的作用机制。

下载地址：http://www.ciejournal.org/Magazine/show/?id=51567

【Data & Stata Files】，右击另存即可。

作者将企业的生命周期划分为三个阶段，分别为初创期、成长成熟期和衰退期, 使用 Cox PH 模型分析了工资扭曲对企业从一个阶段跨越到另一个阶段的影响，并用 AFT 模型做相关稳定性检验。下面截取文章基本回归结果及作者给出的相关代码。

Source: 吴先明等 (2017)，表 3. Data & Stata Files

从初创期到成长成熟期, 工资扭曲 (Wdist) 的系数显著为负, 表明工资扭曲在整体上促进初创企业的成长, 从成长成熟期到衰退期，工资扭曲 (Wdist) 的系数为负但不显著。

对应的 Stata 估计命令如下：

*-吴先明等 (2017) 表 3 - Stata 命令 - 基于 Cox PH 的回归分析
*-完整文件地址：http://www.ciejournal.org/Magazine/show/?id=51567

/////以企业进入成长成熟期作为兴趣事件构建生存分析数据/////
stset year, id(panelid2) failure(life==1) origin(xkysj) 

/////表3基于Cox PH模型的回归分析-从初创期到成长成熟期部分/////
stcox wdist  pdens_100  size  lever   cdens_100  subs roa export state  foreig  ///
 HHI grow  mature  decline , nohr
est store model1

stcox wdist  pdens_100 wdist_pdens_100  size  lever   cdens_100  subs roa export ///
 state  foreig    HHI grow  mature  decline , nohr
est store model2

stcox wdist  pdens_100 wdist_pdens_100  size  lever   cdens_100  subs roa export ///
 state  foreig    HHI grow  mature  decline _Iindu* , nohr
est store model3

stcox wdist  pdens_100 wdist_pdens_100  size  lever   cdens_100  subs roa export ///
state  foreig    HHI grow  mature  decline _I* , nohr
est store model4

stset, clear

/////以企业进入衰退期作为兴趣事件
/////但以企业进入成长成熟期作为进入研究的时间，构建生存分析数据
stset year, id(panelid2) failure(life==2) origin(xkysj) enter(life==1)

/////表3基于Cox PH模型的回归分析-从成长成熟期到衰退期部分/////
stcox wdist  pdens_100  size  lever   cdens_100  subs roa export state  foreig ///
 HHI grow  mature  decline , nohr
est store model5

stcox wdist  pdens_100 wdist_pdens_100  size  lever   cdens_100  subs roa export ///
 state  foreig    HHI grow  mature  decline , nohr
est store model6

stcox wdist  pdens_100 wdist_pdens_100  size  lever   cdens_100  subs roa export ///
 state  foreig    HHI grow  mature  decline _Iindu* , nohr
est store model7

stcox wdist  pdens_100 wdist_pdens_100  size  lever   cdens_100  subs roa export ///
 state  foreig    HHI grow  mature  decline _I* , nohr
est store model8

///导出表3基于Cox PH模型的回归分析
outreg2 [*] using 表3基于Cox PH模型的回归分析.xls, replace

4. 结语

生存分析可以有效地解决数据截堵问题，是研究风险率、生存时间的一种优良的计量方法, 本文通过对生存分析的概念、模型、分析步骤及 Stata 命令等进行介绍，来让大家初步地了解这一方法。考虑到这篇文章是作为生存分析入门使用, 相关模型估计并没有做深入的说明, 这也是本文非常大的不足之处。如果大家想更深入地学习理论, 推荐 Cameron and Trival (2005, Chapter17-19)。

参考文献

Cameron, A.C. and P.K. Trivedi. 2005, Microeconometrics Methods and Applications. Cambriage University Press, New York, NY.
Cleves, M., W. Gould, and R. Gutierrez. 2002. An Introduction to Survival Analysis Using Stata. College Station, TX: Stata Press. [PDF]
Survival Analysis Using Stata - Statistical Horizons. [PDF]
陈强.高级计量经济学及 Stata 应用 [M] . 高等教育出版社, 2014.
陈勇兵,钱意,张相文.中国进口持续时间及其决定因素[J].统计研究,2013,30(02):49-57.
过新伟,胡晓.CEO变更与财务困境恢复——基于ST上市公司"摘帽"的实证研究[J].首都经济贸易大学学报,2012,14(03):47-54.
何树全,张秀霞.中国对美国农产品出口持续时间研究[J].统计研究,2011,28(02):34-38.
李宏伟.基于生存分析的上市公司财务困境时间效应研究[J].财会通讯,2019(32):31-35.
李思瑶,王积田,柳立超.基于生存分析的P2P网络借贷违约风险影响因素研究[J].经济体制改革,2016(06):156-160.
刘海洋,林令涛,黄顺武.地方官员变更与企业兴衰——来自地级市层面的证据[J].中国工业经济,2017(01):62-80.
逯宇铎,于娇,刘海洋.集聚经济是否影响了企业生命周期——基于企业退出行为视角[J].财经科学,2013(10):60-70.
陆志明,何建敏,姜丽莉.基于生存分析模型的企业财务困境预测[J].统计与决策,2007(21):174-176.
吕长江,赵岩.中国上市公司特别处理的生存分析[J].中国会计评论,2004(02):311-338.
毛其淋,盛斌.贸易自由化、企业异质性与出口动态——来自中国微观企业数据的证据[J].管理世界,2013(03):48-65+68+66-67.
毛其淋,许家云.中国对外直接投资促进抑或抑制了企业出口?[J].数量经济技术经济研究,2014,31(09):3-21.
毛其淋,许家云.中国企业对外直接投资是否促进了企业创新[J].世界经济,2014,37(08):98-125.
蒋灵多,陈勇兵.出口企业的产品异质性与出口持续时间[J].世界经济,2015,38(07):3-26.
靳永爱,陈杭,李芷琪.流动与女性生育间隔的关系——基于2017年全国生育状况抽样调查数据的实证分析[J].人口研究,2019,43(06):3-19.
宋雪枫,杨朝军.财务危机预警模型在商业银行信贷风险管理中的应用[J].国际金融研究,2006(05):14-20.
谭晶荣,童晓乐.中国与金砖国家贸易关系持续时间研究[J].国际贸易问题,2014(04):90-100.
王晓鹏,何建敏,马立成.Cox模型在企业财务困境预警中的应用[J].价值工程,2007(11):4-8.
吴冰.生存分析及其应用:以创业研究为例[J].上海交通大学学报(哲学社会科学版),2006(03):63-65+75.
吴先明,张楠,赵奇伟.工资扭曲、种群密度与企业成长:基于企业生命周期的动态分析[J].中国工业经济,2017(10):137-155.
许家云,毛其淋.政府补贴、治理环境与中国企业生存[J].世界经济,2016,39(02):75-99.
张世伟,赵亮.农村劳动力流动的影响因素分析——基于生存分析的视角[J].中国人口·资源与环境,2009,19(04):101-106.
赵瑞丽,孙楚仁,陈勇兵.最低工资与企业出口持续时间[J].世界经济,2016,39(07):97-120.

Appendix 本文涉及的 Stata 代码

use http://www.stata-press.com/data/cggm3/hip, clear //调用 hip.dta 数据
describe //对数据集进行描述
stset time1, origin(time time0) id(id) failure(fracture == 1) //设定数据
stdescribe //在 stset 命令之后使用的, 用于描述生存分析数据的基本特征

*-----使用 sts graph 命令 figure5 -----
use http://www.stata-press.com/data/cggm3/hip2, clear //数据已设置成生存分析数据
sts graph  //生存分析图
graph export figure5.png

*-----使用 sts graph 命令 figure6-----
sts graph, by(protect) plot2opts(lp("-")) //以是否穿防护装置分组画生存分析图  plot2opts(lp("-"))表示第二个图用虚线表示
graph export figure6.png

*-----使用 sts graph 命令 figure7-----
sts graph, cumhaz ci //cumhaz 选项表示画累积风险函数, ci 显示 95% 的置信区间。
graph export figure7.png

*-----使用 sts graph 命令 figure8-----
sts graph, cumhaz ci //cumhaz by(protect)选项表示画累积风险函数, ci 显示 95% 的置信区间。
graph export figure7.png

*-----使用 streg 命令 table5-----
streg protect age calcium, nohr nolog dist(weib)
* dist(weib)代表威布尔回归的风险函数, 
* dist(e)指数回归，
* dist(gom)刚珀茨回归，
* nohr显示回归系数而不显示风险比率

*-----使用 stcurve 命令 figure 9-----
stcurve, hazard //画风险函数图

*-----使用 stcurve 命令 figure 10-----
stcurve, hazard //画风险函数图

*-----使用 streg 命令 table6-----
streg protect age calcium,  nolog dist(logn)

*-----使用 streg 命令 table 7-----
. stcox protect age calcium,r nohr nolog
 
*-----使用 stcurve 命令 figure 11-----
stcurve, hazard //画风险函数图

*-----使用 stphplot 命令 figure 12-----
. stphplot, by(protect)
 
*-----使用 stcoxkm 命令 figure 13-----
. stcoxkm, by(protect)
 
*-----使用 estat phtest命令 table8 -----  
quitely stcox protect age calcium,r nohr nolog
estat phtest, detail

*-----使用 estat phtest命令 figure14 -----   
estat phtest, plot(protect)

专题	嘉宾	直播/回看视频
⭐ 最新专题 ⭐		DSGE, 因果推断, 空间计量等
⭕ Stata数据清洗	游万海	直播, 2 小时，已上线
研究设计	连玉君	我的特斯拉-实证研究设计，-幻灯片-
面板模型	连玉君	动态面板模型，-幻灯片-
面板模型	连玉君	直击面板数据模型 [免费公开课，2小时]