生存分析（3）

忘是亡心i 2022-06-08 03:38 1877阅读 0赞

之前的一篇文章，我们比较细致的讲解了一下K-M法与寿命表法这两种**非参数**的研究方法。实现了对数据的一个初步应用，并能够完成**单因素**分析，接下来将介绍一下生存分析中最常用的模型—Cox比例风险回归模型，这是一种**含参数**的模型，能够实现对数据的**多因素分析**。

#### Cox比例风险回归模型 ####

Cox比例风险回归模型（Cox’s proportional hazards regression model），简称Cox回归模型。  
该模型由英国统计学家D.R.Cox于1972年提出，主要用于肿瘤和其它慢性病的预后分析，也可用于队列研究的病因探索。其优点：  
(1) 多因素分析方法  
(2) 不考虑生存时间分布  
(3) 利用截尾数据

**Cox回归模型是生存分析最重要也最常用的模型，是一种半参数模型**

##### 1.Cox模型的基本形式 #####

h(t,X)=h0(t)exp(β1X1\+β2X2\+⋯βpXp)

h(t,X)是 **风险函数**  
即生存时间已到达t的一群观察对象在t时刻的瞬时死亡率。（ [详见生存分析（1）][1]）

h0(t)是**基准风险率**，是与时间有关的任意函数，函数形式无任何限定（此时所有协变量取值为0）

Xi为**协变量**，即影响因素  
函数右侧分为两部分一个是无规定的h0(t)属于非参部分，另一部分是含参的部分，因此Cox模型又称为**半参数模型**。

βj是**回归系数**  
βj>0时，协变量的取值越大，风险函数h(t,X)的值越大，表示病人死亡的风险越大  
βj=0时，表示协变量对风险函数h(t,X)没有影响  
βj<0时，协变量的取值越大，风险函数h(t,X)的值越大，表示病人死亡的风险越小

某危险因素Xi在非暴露组取值为0，在暴露组取值为1

h(t,Xi=1)h(t,Xi=0)=h0(t)exp(βi)h0(t)=exp(βi)=RR

式中暴露组与费暴露组的危险率之比，称为 **相对危险度RR**  
易知,  
RR>1，表示风险越大  
RR=1，变量对生存无影响  
RR<1，表示风险越小

利用h(t,X)与协变量的关系还可以推导出S(t,x)与协变量的关系

h(t)=f(t)S(t)=F′(t)S(t)=(1−S(t))′S(t)=−S(t)′S(t)

两边同时对t求积分，得  −∫t0h(t)dt=ln(S(t)) ，结果带入回归方程：

S(t,X)=exp\[−∫t0h(t,X)dt\]=exp\[−∫t0h0(t)exp(β′X)dt\]=\[S0(t)\]exp(β′X)

其中  S0(t)=exp\[−∫t0h0(t)dt\]

##### 2.Cox模型的前提条件 #####

Cox模型必须满足比例风险假定（**PH假定**）

*  任何两个个体的风险函数之比，即风险比（HR）保持一个恒定的比例，与时间t无关
 *  模型中协变量的效应不随时间改变而改变

其中，

HR=hi(t)hj(t)=exp\[β1(Xi1−Xj1)\+β2(Xi2−Xj2)\+⋯\+βp(Xip−Xjp)\]

检查某协变量是否满足PH假定，最简单的方法是观察按该变量分组的生存曲线，若**生存曲线交叉**，提示不满足PH假定

PH检验的思想是这样的，令

h(t,X)=h0(t)exp\[(β1\+G(t)θ)X1\+β2X2\+⋯\+βpXp\]

其原假设为H0：  θ=0 ，如果验证了则说明该协变量的效应（系数）与时间无明显关系，所以P值应越大越好。

**不满足PH假定时，可采用分层比例风险模型**

##### 3.参数估计——偏似然估计 #####

**危险集：**对于每个生存时间ti来说，凡生存时间**大于等于ti**的所有病人组成一个危险集。

下面以一个例子说明如何对βj进行偏似然估计  
![这里写图片描述][SouthEast]  
**条件死亡概率**

qi=(ti,ti\+1)时刻死亡的例数ti时刻尚存活的例数=(ti,ti\+1)时刻死亡的例数≥ti时刻死亡或右删失的例数=ti时刻尚存活例数×hi(t)∑j≥i\[ti时刻尚存活例数×hj(t)\]

=hi(t)∑j≥ihj(t)  
因此对应上表中的三个死亡时刻（18、48、90），有三个条件死亡概率如下：

q1=eb1\+b2eb1\+b2\+e0\+eb2\+eb1，q2=11\+eb2\+eb1，q3=eb1eb1

将上述各个条件概率相乘可以得到偏似然函数  Lp=q1q2…qk  
各个条件死亡概率可以表示为：

qi=hi(t)∑j≥ihj(t)=h0(t)exp(β1Xi1\+β2Xi2\+⋯\+βpXip)∑S∈R(ti)h0(t)exp(β1Xs1\+β2Xs2\+⋯\+βpXsp)

=exp(β1Xi1\+β2Xi2\+⋯\+βpXip)∑S∈R(ti)exp(β1Xs1\+β2Xs2\+⋯\+βpXsp)

其中，S表示ti时刻以后危险集R(ti)中对似然函数做贡献的个体  
将n个病人的死亡条件概率相乘得：

L(β)=∏i=1nexp(β1Xi1\+β2Xi2\+⋯\+βpXip)∑S∈R(ti)exp(β1Xs1\+β2Xs2\+⋯\+βpXsp)

第  i 个研究对象在  ti 时刻死亡的概率应当是两部分的乘积，一是患者存活到  ti 时刻的概率（与  h0(t) 有关），二是该暴露人群  Ri 中恰好第i个患者死亡的概率（  qi ），  L(β) 忽略了前者，故称之为 **偏似然函数**

为了体现删失数据不参与计算的特点，用δi=1表示病人在ti时刻死亡，用δi=0表示病人在ti时刻删失，其偏似然函数为：

L(β)=∏i=1n\[exp(β1Xi1\+β2Xi2\+⋯\+βpXip)∑S∈R(ti)exp(β1Xs1\+β2Xs2\+⋯\+βpXsp)\]δi

两边同时以e为底取对数，得：

lnL(β)=∑i=1n\{ δi\[(β1Xi1\+⋯\+βpXip)−ln∑S∈R(ti)exp(β1Xs1\+⋯\+βpXsp)\]\}

求关于  βj（j=1,2,…p） 的一阶偏导数，并使其等于0，即令  ∂lnL(β)∂βj=0  
可解得  βj的最大似然估计  
现在我们已经拥有了计算  βj 的方法，那么接下来要解决的问题是：究竟选择哪些协变量参与模型呢？

##### 4.影响因素的筛选 #####

*  筛选方法有前进法、后退法和逐步回归法、最优回归子集法，实际工作根据具体情况选择使用，最常用的为**逐步回归法**。
 *  因素筛选时需规定显著性水平，一般确定为0.1或0.15，设计较严格的研究可定为0.05。
 *  筛选因素时，还要考虑因素间是否有共线性。若存在共线性，应考虑**消除共线**的影响，如采用主成分回归等方法。

**参数假设检验**  
Cox提供三种基于极大似然法的大样本检验方法，通过检验结果选择协变量。

*  **似然比检验(likelihood ratio test)**
 *  **Wald检验**
 *  **得分检验(score test)（又称为拉格朗日乘数法）**

这三种假设检验的原假设或者约束都是H0:βj=0  
若原假设成立则不选择该协变量，否则选择该协变量进入回归模型

1）似然比检验（LR）  
用于模型中原有不显著变量的剔除和显著变量的引入，以及包含不同变量数时模型间的比较。  
检验新增加的协变量是否有统计学意义的统计量为：

χ2=2ln\[lnL(p\+1)−lnL(p)\]

χ2 服从自由度1的  χ2 分布，如果差值无统计学意义则不引入改协变量，反之则引入。  
把原有模型中无统计学意义的协变量剔除，其方法与增加协变量的检验方法相似。

2）Wald检验  
它用于模型中的协变量是否应从模型中剔除。检验模型中协变量Xj对模型的贡献是否有统计学意义，其对应的wald统计量为:

χ2ω=(βjSβj)2

它服从自由度为1的  χ2 分布。式中  Sβj 表示回归系数  βj 的标准误  
当  βj 的95%可信区间包含0时，  βj 为0

3）拉格朗日乘子检验（LM）

*  用于一个或多个新变量是否能入选模型
 *  用于检验变量间的交互作用能否对生存时间产生影响

**基本思想：**拉格朗日乘子检(LM)，又称为Score检验。该检验基于约束模型，无需估计无约束模型。

假设约束条件为H0:g(θ)=C，在约束条件下最大化对数似然函数，令λ表示拉格朗日乘子向量，此时，拉格朗日函数为：

LnL∗(θ)=LnL(θ)\+λ\[g(θ−C)\]

约束条件下最大化问题就是求解下式根，  
 ∂Ln∗L(θ)∂θ=∂LnL(θ)∂θ\+λg′=0  
 ∂Ln∗L(θ)∂λ=g(θ)−C=0

如果约束条件成立，对数似然函数值不会有明显变化。这就意味着在一阶条件下，第二项应该很小，特别是λ应该很小。因此约束条件是否成立检验转化成检验H0:λ=0，这就是拉格朗日乘子检验的思想.

但是直接检验H0:λ=0比较困难，有一个等价而简单的方法。  
如果约束条件成立，在约束估计值处计算对数似然函数的导数应该近似为零，如果该值显著异于零，则约束条件不成立，拒绝原假设

对数似然函数的导数就是得分向量，因此，LM检验就是检验约束条件下参数估计值的得分向量值是否显著异于零，因而，LM检验又称为得分检验。

**三种假设检验**  
LR:需要同时最小化有约束和无约束模型  
Wald:只需要优化无约束模型  
LM:只需要约束有约束模型

LR、Wald、LM关系（一般情况下成立）：

Wald≥LR≥LM

##### 5.Cox模型+Roc预测生死？ #####

这是基于我所看到的一些医学文献的做法做的小小总结，实话说小生之前对ROC不是很明白，为此看了不少内容。之后应该会在统计类博客里面做更详细的介绍，只想了解大概的同学可以看[这儿][Link 1]

一般来说通过Cox模型我们能够得到风险函数，也就能够通过每个个体的特征计算其“风险值”（与存活时间t相关）。但风险值是一个回归结果，只能表示患者死亡的风险大小并不能直接用来做分类，我们必须找到一个“风险值”的临界点，当患者的“风险值”大于这个临界值则预估其在t时刻死亡，否则存活。

一般来说ROC曲线是用来描述二分类模型的优劣的，而为了描绘曲线需要多个点，如何通过一个Cox模型的风险函数获得多组分类结果呢？

我们可以通过将“风险值”从小到大排序，将每一个“风险值”都作为一次临界值，分别以此判断患者的生存情况并计算AUC，AUC最大时的“风险值”则是最优临界值，我们便可以以此为临界值去预估其他患者在t时刻的生存情况啦~~~

**算例：**  
我们根据每个测试样本属于正样本的概率值从大到小排序。下图是一个示例，图中共有20个测试样本，“Class”一栏表示每个测试样本真正的标签（p表示正样本，n表示负样本），“Score”表示每个测试样本属于正样本的概率。

![这里写图片描述][SouthEast 1]

接下来，我们从高到低，依次将“Score”值作为阈值，当测试样本属于正样本的概率大于或等于这个阈值时，我们认为它为正样本，否则为负样本。举例来说，对于图中的第4个样本，其“Score”值为0.6，那么样本1，2，3，4都被认为是正样本，因为它们的“Score”值都大于等于0.6，而其他样本则都认为是负样本。每次选取一个不同的阈值，我们就可以得到一组FPR和TPR，即ROC曲线上的一点。这样一来，我们一共得到了20组FPR和TPR的值，将它们画在ROC曲线的结果如下图：

![这里写图片描述][SouthEast 2]

灵敏度-特异度最大时的“阈值”便是我们的最优临界值。

##### 6.几种参数回归模型 #####

生存分析中常见的回归参数模型，主要有**指数**回归模型、**Weibull**回归模型、**对数logistic**回归模型、**Gamma**模型等。

**1.指数回归模型**  
指数分布是历史上第一个生存时间分布模型,是最简单与最重要的分布,在生存分析中占据重要地位。指数分布是一种**纯随机死亡模型**,在任何时间上的风险函数为一常数,即风险函数的大小不受生存时间长短的影响,以独特的“无记忆性”而闻名,常被看成随机失效(死亡)模型。  
设数据来自指数分布，其**概率密度函数**为：

f(t)=\{ λe−λt,0,t≥0,λ>0t<0

**分布函数**为：

F(t)=1−e−λt (t≥ 0)

**生存函数**为：

S(t)=1−F(t)=e−λt

**风险函数**为：

h(t)=f(t)/S(t)=λλ>0,t≥0

λ 为指数分布的风险率,称为刻度参数或尺度参数,其大小决定了生存时间的长短。风险率越大,生存率下降越快,生存时间越短;风险率越小,生存时间越长。在指数分布模型中,  λ 为常数,与时间T无关,这是指数分布的特点,可利用此特点识别某生存时间分布是否服从指数分布。  
指数回归是在指数分布的基础上引入预后因素后的模型,指数回归模型在疾病的影响因素研究中起很重要的作用,在疾病的影响因素研究中,患者的生存状况总是和其他的某些指标以及病情的某些特征紧密相关的,因此,在考虑某一因素对病人生存时间的影响时,还必须考虑其他影响因素或协变量的影响,并对其作用大小进行分析,需建立多因素回归模型,从而查明各种因素如何影响生存率。  
在指数回归模型中,设  x1,x2,…xp 为 影响因素，如果生存时间服从指数分布，则参数  、lambda 与各因素间的关系可用如下回归方程表示：

ln(λ)=β0\+β1x1\+β2x2\+⋯\+βpxp

式中  β0 为常数项，表示无任何因素影响时的基准风险之对数。（比一般Cox模型多了一个  β0 ）  βi 表示控制其他影响因素时，变量  xi 每改变一个单位所引起的风险函数  λ 之对数值的该变量。风险函数为：

h(t)=λ=exp(β0\+β1x1\+β2x2\+⋯\+βpxp)

基准风险表示其余因素都“不存在”的情况下的风险。

h0(t)=λ0=eβ0

相应t时刻的生存率为：

S(t)=exp\[−t∗exp(β0\+β1x1\+β2x2\+⋯\+βpxp)\]

**2.Weibull模型**  
Weibull分布也是生存分析的理论基础,由瑞典科学家Waloddi Weibull提出。Weibull分布是一种连续性分布，它是指数分布的一种推广形式,它不像指数分布假定危险率是常数,因而有更广的应用性。  
f(t)的**概率密度函数**为：

f(t)=λγ(t)γ−1exp\[−λ(t)γ\]

**分布函数**：

F(t)=1−exp\[−λ(t)γ\]

**生存函数**：

S(t)=exp\[−λ(t)γ\]

**风险函数**：

h(t)=λγ(t)γ−1

其中  λ 和  γ 为两个参数。  λ 称为 **尺度参数**,它决定分布的分散度;  γ 为 **形状参数**,它决定该分布的形态。  γ>1 时风险函数随时间单调递增;  γ<1 时风险函数随时间单调递减;显然,当  γ=1 时,风险不随时间变化,weibull分布简化为指数分布,所以指数分布是weibull分布在  γ=1 时的特例。  
如果资料服从Weibull分布,则可用回归模型对危险因素进行分析。在Weibull回归模型中,风险函数与影响因素间的关系也假设为指数关系,即:

λ=exp(β0\+β1x1\+β2x2\+⋯\+βpxp)

**风险函数**为：

h(t)=γ(t)γ−1exp(β0\+β1x1\+β2x2\+⋯\+βpxp)

**生存函数**为：

S(t)=exp\[−tγexp(β0\+β1x1\+β2x2\+⋯\+βpxp)\]

**基准风险**为：

h0(t)=γ(t)γ−1eβ0

**3.Gamma模型**  
生存分析讨论两类不同的Gamma模型:标准Gamma模型(2参数)和广义Gamma模型(3参数)。  
标准Gamma分布的特性取决于两个参数γ和λ,γ为形状参数,λ为尺度参数。当0<γ<1时,若时间从0增加到无穷时,风险函数从无穷单调地减小到λ,表现为负老化;当γ>1时,若时间从O增加到无穷时,风险函数从0增加到λ，表现为正老化；当γ=1时,风险等于常数λ，即指数分布情形。  
广义Gamma模型可以呈现不同于其他任何分布特例的形状。特别地，它可以是U形或浴盆形的风险函数，在这样的函数中风险函数先下降，下降到最低值后又升高。总所周知，人累在整个生命周期中的死亡危险性就属于这种形状。

有关实用软件实现生存分析，由于编者比较熟悉R，因此只介绍R的实现。具体内容，请参考R语言学习类博文。

[1]: http://blog.csdn.net/xiaohukun/article/details/77679134
[SouthEast]: /images/20220608/aa440aef35f2491594e4cd8554bab23b.png
[Link 1]: https://www.douban.com/note/284051363/
[SouthEast 1]: /images/20220608/3637d9b490e542449ebe76255f5db8b4.png
[SouthEast 2]: /images/20220608/192a990187ce40edb9df7cc7d1a6046b.png