生存分析-基础概述和常用方法

最新推荐文章于 2024-07-09 23:41:00 发布

No..D

最新推荐文章于 2024-07-09 23:41:00 发布

阅读量1.3k

点赞数 25

分类专栏：笔记文章标签： r语言数据分析

本文链接：https://blog.csdn.net/weixin_45275365/article/details/138589157

版权

笔记专栏收录该内容

8 篇文章 1 订阅

订阅专栏

碎碎念

我导有和医院合作的项目意向，提前补充一点医学上做生存分析的基础概念，相关方法和R语言实现知识。

1. 基础知识

1.1 基本概念

生存：广义上指的是研究对象是否出现感兴趣的阳性终点事件，或者说是某一现象是否出现失效事件（failure），例如研究复发情况时，复发=死亡，未复发=生存。

生存数据：time-to-event，因变量由两个部分刻画，结局状态与生存时间，即终点事件是否发生以及发生终点事件所经历的时间。

生存分析：包括生存时间的分布描述、生存时间分布的组间比较、生存时间分布的影响因子以及效果评估。核心目标是研究群体的死亡速度，又或者是规定时间内的死亡率或生存率。

1.2 数学抽象

生存概率 $S (t)$ ：研究对象从开始到某个时间点t仍然存活的概率

$S(t)=P(T\geq t)$

随着时间增加，S(t)减小，因此$S(0)=1,S(\infty)=0 $

风险概率 $h (t)$ ：研究对象从开始到某个时间t之前存活，但是在t时间点死亡的概率

$h(t)=-\frac{d\ln S(t)}{d t}=\lim_{\Delta t \rightarrow 0}{\frac{P(t\leq T < t+\Delta t | T\geq t)}{\Delta t}}$

累积风险函数 $H (t)$ ：也就是 $h (t)$ 从0开始到一定时间内的积分

$H(t)=\int_{0}^{t} h(u)du=-\ln S(t)$

1.3 研究目的分类

生存时间的分布描述：根据收集到的生存资料估计总体生存率或者其他的有关指标。例如：获得某种肿瘤患者不同时间点的生存率、总体的生存曲线以及中位生存时间等等。

生存时间分布的组间比较：对于不同处理组，分析不同组之间的生存率、生存时间等指标。例如：根据不同的治疗方案将患者分组，比较不同组之间的生存曲线，以了解哪种方案较优。

生存时间分布的影响因子：探索影响生存时间长短的相关因素，又或者固定某些因素后，研究其他因素对生存率的影响。例如：探索影响癌症患者预后的主要因素，包括病人的性别、肿瘤分期、治疗方案等等。

效果评估：根据不同因素水平预测患者个体生存概率。

1.4 常见临床终点

OS，Overall Survival：总生存期，治疗结束 → 患者死亡

DFS，disease-free survival：无病生存期，治疗结束 → 肿瘤复发或者任意原因死亡

PFS，progression-free survival：无进展生存期，治疗结束 → 肿瘤客观进展或者患者死亡

DMFS，distant metastasis-free survival：无远处转移生存期，治疗结束 → 发现有转移病灶

LRRFS，locoregional relapse-free survival：局部区域无复发生存期，治疗结束 → 发现局部-区域复发

2. 常用方法

2.1 Kaplan-Meier，KM曲线

一种非参数方法，用于生存率的估计和影响因素分析
图片来自于文献1 思路：假设研究时间点t时刻的S(t)，即计算一个病例在t时刻发生死亡的概率P(T≥t)，那么逻辑上，在t之前的每个观察时刻，该病例都应该存活，即发生（1-P(T≥t’)）的事件，直到t时刻发生死亡。因此根据每个观察时刻的情况，计算死亡概率，并利用累乘计算，即可获得累积生存概率函数S(t)。

数学表示：m个观察时间点，生存概率可表示为

$S(t_i)=S(t_{i-1})(1-\frac{d_i}{n_i})\\ d_i为i时刻的死亡人数，n_i为i时刻的有效存活人数$

比较分析：利用统计学工具量化指标衡量不同KM曲线之间的差异和显著性。

中位生存时间：计算生存率为50%时所对应的时间。
Logrank对数秩检验：一种非参数检验方法（卡方检验），比较两组样本生存时间分布的差异。

零假设：两组生存时间分布完全一致

计算某个时间点期望死亡人数 $E_{it}=N_{it}\cdot (O_t / N_t)$ ，O为实际死亡人数，N为有效风险人数，i表示为第i组对应指标。如果零假设成立，则后面的人数占比，全局和第i组应该相同。

构造卡方值： $\chi ^2=\sum_j \frac{(\sum O_{jt}-\sum E_{jt})^2}{\sum E_{jt}}$ ，然后利用查表确定是否拒绝零假设
Breslow假设检验：非参数检验方法（Wilcoxon检验），计算不同时间点上的期望人数时有一个权重因子（通常为当前的有效风险人数），然后再去构造卡方值进行检验。

置信区间：每个观察点相当于做了一次随机抽样，这意味着存在误差，因此可以利用S(t)获得特定概率的置信区间，以确定误差。

$S(t)\pm z_{\alpha /2}\sqrt{Var(S(t))}\\ Var(S(t))=S^2(t)\sum\frac{d_i}{n_i(n_i-d_i)}$

该公式的问题是，假设存在ni=di的情况，会导致求和的最后一项失效。因此有新的公式：

$\exp(-\exp(c_{\mp}))\\ c_{\mp} =\log(-\log S(t)) \mp z_{\alpha/2}\sqrt{V}\\ V=\frac{1}{(\log S(t))^2} \sum \frac{d_i}{n_i(n_i-d_i)}$

2.2 Logistics分析

因变量只考虑分类，也就是事件发生情况，根据线性回归模型发展而来，一种参数类模型。

数学公式：

$P(y=1)=\frac{\exp(\beta_0+\sum_i \beta_ix_i)}{1 +\exp(\beta_0+\sum_i \beta_ix_i)}=(1+\exp(S))^{-1}$

其中$\beta $为偏回归系数，表示在其他自变量固定的情况下，某个自变量每改变一个单位，log§的平均改变量。该数值与OR值相关：

$OR_j=\exp[\beta_j(c_1-c_0)]$

OR指的是比值比，表示某个因素与事件发生的相关性，OR=1，表示不相关，OR>1，表示正相关，OR<1，表示负相关。例子： $OR=\frac{肺癌中吸烟/肺癌中不吸烟}{健康中吸烟/健康中不吸烟}$

参数估计：偏回归系数常常利用极大似然估计方法进行估计。

检验假设：

模型检验：零假设为 $\beta_0=\beta_1=\beta_2=...=\beta_m=0$ ，采用似然比G检验等方法进行。
参数检验：零假设为 $\beta_i=0$ ，通常采用卡方检验。
拟合优度检验：即检验模型获得数据与实际的符合情况

Cox分析

比例风险模型，一种半参数的分析方法。针对time-to-event数据。不直接对S(t)建模，而是对风险概率h(t)建模。

数学公式：

$h(t)=h_0(t)\cdot \exp(\sum_j \beta_jx_j)\\ h_0为基准风险函数$

基准风险函数通常可以利用非参数模型确定，而后续的部分则是利用参数β确定。

参数估计：采用极大似然函数确定。

参数解析：

通常会获得某个协变量的不同风险比（hazard ratio，HR） $HR=\exp(\beta)$ ，HR=1意味着这个因素与结果之间不存在作用。HR>1证明该变量增加，则风险增加，HR<1则说明该变量减小，风险增加。

前提假设：任意两人的风险比是不随着时间变化而变化的。

对两组数据做KM曲线，曲线如果交叉，则说明两组的生存概率关系随着时间变化而发生了变化。
以ln(t)为横坐标，ln(-ln(S(t)))为纵坐标，两组之间应该近似平行或者等距离
Schoenfeld残差检验

2.3 倾向得分匹配（Propensity score matching）

为了研究某个因素对结果的影响，通常情况下，会利用该变量进行分组，对比组别之间的差异。举例来说，研究吸烟与肺癌的复发情况，收集了一队列数据，对比分析吸烟人群与不吸烟人群的复发差异。但是我们实际上是想知道，对于病人A，他吸烟和不吸烟，复发的概率差异有多少。然而病人A吸烟已经是事实，那怎么样可以估计出病人A不吸烟的时候的复发概率呢，这就是倾向得分匹配面对的问题。

基本步骤：

建立倾向性评分：倾向性评分指的是，已知协变量取值为 $X_i$ 的个体，受到干预的概率 $E (X) = P (G = 1∣ X)$ 。假定个体所在组别与协变量无关，即G与X相互独立。可以用logistic回归或其他方法计算。
匹配受试者：接受干预与未干预的个体之间，寻找倾向得分相似的匹配对。
评估匹配效果：确保匹配后的两组在干预特征上更平衡，也就是保证匹配后两组在关键特征上更加相似
比较结果：匹配后，比较干预组与对照组的结果。

也就是大致可以划分为计算倾向得分和利用匹配策略匹配两部分。

R语言实现

library("tableone")
library("MatchIt")
library("survival")
library("survminer")

# 读取数据
dat <- read.spss('D:\project\0_other\learnR\seer白血病\seercoxleukaemia.sav', to.data.frame=T)
# 数据转换
dat$age_grp <- ifelse(dat$age<=60,1,2)
dat$age_grp <- factor(dat$age_grp, labels=c('<=60', '>60'))
# 批量进行因子化
vars<-c("Sex","Race","matrimony","Chemotherapy",
        "Radiation","Grade","Stage","age_grp","css2")
dat[vars]<-lapply(dat[vars],factor)

# 显示基础的数据结果
myVars <- c("age","Sex","Race","matrimony","Chemotherapy",
            "Radiation","Grade","Stage","status1",
            "css2","age_grp","Year.of.diagnosis")
catVars <- c("Sex","Race","matrimony","Chemotherapy",
             "Radiation","Grade","Stage","status1",
             "css2","age_grp")
tab1 <- CreateTableOne(vars = myVars, data = dat, factorVars = catVars,
                       strata = "Chemotherapy")
# vars: 需要描述的统计量，strata: 分组变量，data: 设定数据集，factorVars: 指定分类变量
print(tab1) # 匹配前的数据

# PSM
formula1<-Chemotherapy ~ age+Sex+Race+matrimony+Radiation+Grade+Stage+Year.of.diagnosis
m_out<-matchit(formula1, data = dat, method = "nearest",
               distance="glm", m.order= "random",ratio=1,
               caliper=0.2,replace=FALSE)
# method: 匹配方法，distance：获得评分的方法，ratio: 允许以1：n的比例进行匹配，caliper: 卡钳值，越小越严格

# 评估
psm_matchit_data <- match.data(m_out) # 生成匹配后的数据
tab2 <- CreateTableOne(vars = myVars, data = psm_matchit_data, factorVars = catVars,
                       strata = "Chemotherapy")
print(tab2)

3. 结果分析方法和实现

3.1 Nomogram 列线图

将多因素分析结果，用图形的方式表现出来。

思路：简而言之就是根据模型中各个因素对结局的贡献程度，给每个影响因素的不同取值水平赋分，最后相加得到总的评分，在进行转换。

R语言实现

Logistic回归模型

nomogram参数说明

# 加载所用的库
rm(list = ls())
library(survival)
library(survminer)
library(rms)

## 加载数据
lung$sex <- factor(lung$sex,
                   levels = c(1,2),
                   labels = c("male", "female"))
## 打包数据
dd <- datadist(lung)
options(datadist='dd')

## 构建模型
log_fit <- lrm(status ~ age+sex, data=lung)

## 画图
nom <- nomogram(log_fit, fun=function(x)1/(1+exp(-x)), lp=F, funlabel='Risk')
plot(nom)

log_fit：即构建的模型
fun：用于变换预测值，也就是从总体评分转换为所需要的值的函数，以目前为例，log_fit的输出就是x，要将它转换为风险概率，则需要利用log变化。

在这里插入图片描述
Cox回归模型

## Cox
cox_fit <- psm(Surv(time, status) ~ age+sex, data=lung, dist='lognormal')
med <- Quantile(cox_fit) # 计算中位生存时间
surv <- Survival(cox_fit) # 构建生存概率函数

## lung数据中time以天为单位
nom_cox <- nomogram(cox_fit, fun=list(function(x) surv(365, x),
                                      function(x) surv(730, x)),
                    funlabel = c('1-year Survival Probability',
                                 '2-year Survival Probability'))
plot(nom_cox)

在这里插入图片描述

3.2 Calibration curve 校准曲线

是对模型评估的基础指标之一，用于评价一致性/校准度，即预测值与真实值之间的差异程度。实质上是Hosmer-Lemeshow拟合优度检验结果的可视化。

基本思路：

根据预测模型计算个体未来发生结局事件的预测概率
根据概率值从小到大排序，按照十分位等分为10组
分别计算各组的实际观测数和模型预测数， $N_{模型预测}=\sum_{i} p_{i}*N_{i}$
根据每组实际观测数和模型预测数计算卡方值（自由度=8），在计算对应的P值

P＜0.05，表示预测与实际有一定的差异，代表模型的一致性较差；P>0.05表示两者之间没有显著的差异，模型的一致性较好。

绘制步骤：

校准曲线就是绘制实际观测数和模型预测数的散点图。步骤为：

对预测概率进行分桶（策略：‘uniform’, 'quantile’分位点）
计算每个桶里所有样本预测概率平均值，作为横坐标
求每个桶里正例的概率，作为纵坐标
连线

R语言实现：

cox_fit <- cph(Surv(time, status) ~ age+sex, data=lung, surv=T, x=T, y=T, time.inc=365)
# 绘制校准曲线的拟合需要提前设定x=T, y=T
cal <- calibrate(cox_fit, cmethod='KM', method='boot', u=365, m=38, B=228)
# u=time.inc, m为每组样本数量，3*m大约等于总样本量，B为最大再抽样样本量
plot(cal,
     lwd=2, # 线宽
     lty=1, # 线型
     conf.int=T, # 是否显示置信区间
     errbar.col='blue', # 误差棒颜色
     col='red',# 曲线颜色
     subtitles = F,#副标题是否显示
     xlab = 'Cox-Predicted Probability of 1-year DFS',
     ylab = 'Actual 1-year DFS'
     )

在这里插入图片描述 C-index：Concordance index，一致性指数，用于计算预测值与真实值之间的区分度（Cox分析）。其思路是将所有研究对象随机两两配对，如果真实生存时间长的个体预测的生成时间也长，或者生存概率高的预测生存概率也高，则称之为一致。C-index就是统计这些一致的结果占总配对的比例。C-index在0.5-1之间。

3.3 Forest Plot 森林图

Meta分析的一种图形表示方法，直观显示各项之间的关联。延伸一下，可以用于Cox分析后显示HR的关系。在Cox分析后，每个变量可以用一条横线表示，然后会有一个HR和对应的p值，然后这些信息可以具象化表示为森林图。

图像说明：
图像来自于文献1

竖线表示HR=1，表明研究因素与结局无关联。
方块为每一项的HR点估计值。
横线代表HR的95%置信区间。与灰色线的位置不同表达了不同的信息：
- 交叉：无统计学意义
- 左侧：负向相关，即变量数值越高，风险越低
- 右侧：正向相关，即变量数值越高，风险越高

R语言实现

cox_fit <- coxph(Surv(time, status) ~ age+sex+ph.ecog, data=lung)
ggforest(cox_fit, data=lung,
         main = 'Hazard ratio of Lungs', # 标题
         cpositions = c(0.05, 0.15, 0.35), # 列距离
         fontsize = 1, # 字体大小
         refLabel = 'Reference', # 分类变量的参考标签
         noDigits = 3 # 小数点位数
         )

在这里插入图片描述

3.4 Receiver Operating Characteristic ROC曲线

非常常用的预测模型性能评价工具。根据X对Y的预测值与真实值之间的预测准确率情况进行绘图。

基础概念

混淆矩阵
真阳性率TPR：所有阳性群体被检验出来的概率，即 $TPR = TP / (TP + FN)$ ，越接近1越好，也就是常说的Sen（敏感性）
假阳性率FPR：所有阴性样本被检验为阳性的概率，即 $FPR = FP / (FP + TN)$ ，越接近0越好，值等于1-Spe（特异性）

总的来说，一个好的模型应该是TP与TN越多越好，呈现出来的结果就是TPR→1，FPR→0

思路：ROC曲线，就是以TPR为纵轴，FPR为横轴绘制的图像。通常情况下模型会给每个样本一个得分（理解为预测为阳性的概率），那可以在得分范围内，取若干个得分阈值（cut-off value），每个阈值对应一个混淆矩阵，也就获得一组(TPR, FPR)。将所有的观察点连接起来，就能获得ROC。ROC越靠近左上角越好（TPR→1，FPR→0）。

相关应用

利用ROC分析可以获得最佳的cut-off value：保证TPR高的时候，FPR小，可以建立模型max(TPR+(1-FPR))，求解并获得最佳阈值。可用约登指数衡量：Youden index = Sen + Spe - 1.
AUC值：Area Under Curve，定量分析ROC优劣，越大越好。
性能比较-Delong test

time-dependent ROC

在生存分析中，因变量与自变量都有可能随着时间变化而变化，标准ROC无法考虑时间因素，由此衍生出时间依赖性ROC。

常用的定义模式为cumulative/dynamic（C/D），累积动态。对于任意的时间t，可以划分阳性（ $P_t$ ）与阴性组（ $N_t$ ）。病例A在0-t时间段内出现事件，则划分到阳性组。相应的敏感性和特异性指标计算如下：

$Sen(c,t)=P(X_i>c|T_i\leq t)\\ Spe(c,t)=P(X_i \leq c|T_i>t)$

采用R语言包可以随着时间推进，画出一系列的ROC曲线。

3.5 Decision Curve Analysis 决策曲线（DCA）

Decision Curve Analysis，一种评估临床预测模型的新的简单方法。优势在于考虑临床实际效用问题，比较不同模型之间的优劣。

问题说明
（这个部分我学得有些晕乎，欢迎指正！例子来源于文献2）

根据一个临床场景说明，PSA水平上升是早期前列腺癌高风险的常用指标，但是确定是否是高风险癌症还需要进一步活检才能确定。大多数的PSA水平上升患者，要么是无癌症，要么是低风险，不需要进一步治疗。而活检是一项侵入性的，不舒适的有可能引起感染的项目，也就是，如果根据一个人的PSA水平上升，就假定他是高风险癌症人群，再经过活检才能确定其是否真的是，如果不是，那这个误判就导致了损失。假设有新的研究，利用血液检查指标确定高低风险人群，进一步进行后续的活检。

放入一个具体的例子，可以将它抽象化为：100人的队列，其中25例为阳性（癌症高风险），75例为阴性（低或者无风险），现在有模型A（PSA水平升高与否），模型B（血液检查指标）对高低风险进行分类。获得混淆矩阵如下：漏检：FN，误检：FP。

在这里插入图片描述
临床医生可接受的程度是每发现1位患者是真阳性，能接受不超过10名患者预测为阳性。→ 漏检的危害是误检的9倍。所以我们希望重视漏检的危害。定义net benefit，净收益：Benefit - (harm x exchange rate)。随着exchange rate的变化，净收益也会随之变化，甚至出现相对关系变化。

进一步，可以获得不同的exchange rate的净收益，从而获得决策曲线：

确定threshold probability（ $p_t$ ）→ 获得预测标签
计数：记录TP和FP
总样本量为N，计算net Benefit： $Net\ Benifit=\frac{TP}{N}-\frac{FP}{N} \times \frac{p_t}{1-p_t}$
循环前两步，获得数据对 $p_t, NB_t)$