TCGA系列学习笔记（6）多因素生存分析 – sci666

昵称44608199 2022-05-24 发布于浙江

展开全文

TCGA系列学习笔记（6）多因素生存分析

内容目录

前言1. 背景知识1.1 回顾1.2 生存分析种类1.3 Cox回归1.3.1 模型建立1.3.2 基本假定1.3.3 偏回归系数的意义1.3.4 参数估计与假设检验1.3.5 独立危险因素1.3.6 多因素怎么选择因素去做？2. 用R语言进行Cox回归分析2.1 载入数据2.2 单变量计算COX模型2.2.1 性别单因素结果分析2.2.2 批量单因素分析结果2.3 多因素Cox分析2.3.1 多因素分析结果2.4 可视化展示

前言

久违了~

Reference：

一些生存分析相关基础概念：https://www.jianshu.com/p/1a8ee973b45f
关于删失问题详解：https://www./method_topic_article_detail/300/
Kaplan-Meier曲线原理详解：http://www.360doc.com/content/17/0626/11/6175644_666623573.shtml
log-rank检验和Wilcoxon检验的区别：https://mp.weixin.qq.com/s/XpPpOpeNcIDXbd6es5VnvA
生存分析系统讲解：https://www.jianshu.com/p/559d4a966900
用R来做KM生存分析详解：https://www./archives/647/
Cox与KM生存分析及结果解读：https://www./article/1138
强烈推荐——sthda官网R代码：http://www./english/wiki/cox-proportional-hazards-model
Cox回归生存分析 – 简书：https://www.jianshu.com/p/e80eb4168043
R语言实现及结果解读：https://www./article/1132
【8文合集】全面了解单因素分析和多因素分析：https://www./method_topic_article_detail/583/?ty=methods

最后一个深度好文，主要是从统计学角度解释了单因素分析和多因素分析的结果理解，一定要好好看！再放一个链接：

https://www./method_topic_article_detail/583/?ty=methods

要去看哟~

1. 背景知识

1.1 回顾

先简要回顾下TCGA系列学习笔记（5）单因素生存分析中的一些概念：

生存函数：S(t,X) 表示观察对象的生存时间T大于某时刻t的概率，称为生存函数，又称为累积生存率。
死亡函数：F(t,X)=1-S(t,X)，当观察随访到t时刻的累积死亡率。
死亡密度函数：f(t,X)=对F(t,X)求导，某时刻t的瞬时死亡率，称为死亡密度函数。
危险率(风险)函数：h(t,X)=f(t,X)/S(t,X)，某时刻 t 的瞬时死亡速率除以 t 时刻的存活人数（实际上是一个条件瞬间死亡率）。

Kaplan-Meier curves 与 logrank test tests属于单因素分析的例子，他们研究的是单一变量与生存的关系，并且Kaplan-Meier 与log-rank检验只适用于分类变量，却并不适用于数值型变量，比如我们常见的基因表达。

1.2 生存分析种类

生存分析的方法一般可以分为三类：

参数法：知道生存时间的分布模型，然后根据数据来估计模型参数，最后以分布模型来计算生存率。一般不用，因为目前认为我们不知道生存时间符合什么模型。
非参数法：不需要生存时间分布，根据样本统计量来估计生存率，常见方法Kaplan-Meier法（乘积极限法）、寿命法。而对于Kaplan-Meier法来说，其中的p值我们常用log-rank检验和Wilcoxon检验去求。
半参数法：也不需要生存时间的分布，但最终是通过模型来评估影响生存率的因素，最为常见的是Cox回归模型

备注：

单个变量的Cox回归和K-M法结果不一致时，此时我们还是应该选择Cox的结果，因为参数检验效力高于非参数检验。
多个变量的Cox回归中单个变量的显著性与K-M法不一致，此时我们应该选择cox的结果，因为K-M发只考虑单个变量，而Cox考虑了多个变量。

1.3 Cox回归

Cox比例风险回归分析的优势在于可以分析分类变量与数值变量，并且将生存分析的范围由单因素拓展到多因素的分析。

Cox比例风险回归模型可以同时分析几个因素的生存。另外，统计模型提供每个因素的效果大小。

举个例子：假设要比较两组患者，一组基因型是AB和，一组是ab的患者。如果其中AB组都是较老的个体，而ab组都是较年轻的个体。则存活率的任何差异可归因于基因型或年龄或同时两者。因此，在研究与任何一个因素相关的生存时，通常需要调整其他因素的影响。

1.3.1 模型建立

不能直接以S(t,X)为因变量，以X为自变量进行线性回归的原因：

生存分析研究中的数据包含删失数据。
时间变量t通常不满足正态分布和方差齐性。

Cox比例风险回归模型：以风险率函数h(t,X)为因变量，建立指数形式的回归方程（公式1）：

h(t,X) = h₀(t) exp( β₁X₁ + β₂X₂ + β₃X₃ + … + β_mX_m )

β 称为自变量的偏回归系数，用来评价每个危险因素 X_i的效应大小。

h₀(t)称为基准危险率，是当所有因素都不存在时的危险率。

每个exp( β_iX_i ) 称为HR（hazard ratios）：HR>1是危险因素；HR=1无影响；HR<1有利因素。

Cox回归模型的优点：

未对h₀(t)作任何假定，因此Cox回归模型在处理问题时具有较大的灵活性，属于非参数检验。
在许多情况下，即使h₀(t)未知，仍可以估计出参数β。（意义：将时间和危险因素对风险函数的影响分开，仅仅研究危险因素对生存的影响，而不关心时间对生存的影响）

上述公式可以转化为广义线性模型，以便计算和解释（公式2）：

ln[ h(t,X) / h₀(t) ] = lnHR = β₁X₁ + β₂X₂ + β₃X₃ + … + β_mX_m

引入了一个相对危险度HR值来表示 h(t,X) / h₀(t) 的结果。

1.3.2 基本假定

比例风险假定

各危险因素的作用不随时间的变化而变化（公式3）：

HR = h(t,X) / h₀(t)

不随时间的变化而变化（因为时间这个参数已经被消除了）。因此，Cox比例风险回归模型又称为比例风险率模型(PH Model)。这一假定是建立Cox回归模型的前提条件（审稿人很关注的问题，一般要进行检验）。

对数线性假定

对数风险比应与模型中的自变量应与呈线性关系（公式2）：

ln[ h(t,X) / h₀(t) ] = lnHR = β₁X₁ + β₂X₂ + β₃X₃ + … + β_mX_m

这两个假设需要我们进行检验，具体检验的教程后面会再写，需要的可以自己先去学习，地址如下：

http://www./english/wiki/cox-model-assumptions

1.3.3 偏回归系数的意义

偏回归系数β解释了危险因素Xi 对相对危险度RR的影响。

用另外一个通俗的例子来说明：

统计学最后的成绩由期末考试成绩（占比75%）、平时作业成绩（占比20%）和上课签到（占比5%）三部分组成，于是期末考试的成绩在很大程度上决定了你期末最后成绩的高低，而错过了上课签到并不会让你的成绩降低太多。

类比Cox回归，也就是期末考试成绩的“回归系数”比签到的“回归系数”大，更能左右期末考试成绩的高低。是不是很好理解了呢？

1.3.4 参数估计与假设检验

其他自变量不变的情况下，变量 Xi 每增加一个单位，相对危险度RR的(1-α)置信区间为：

Sβ表示β的标准误。

对于回归模型的假设检验通常采用似然比检验、Wald检验和记分检验，其检验统计量均服从卡方分布，其自由度为模型中待检验的自变量个数。注意，这个是对模型整体的检验。

1.3.5 独立危险因素

这个其实又要涉及到比较多的知识点了，首先是——混杂因素：

如果我们想研究某种药物对AML患者的预后影响，那么患者的年龄、性别、有没有基础疾病等都可能是这个研究的混杂因素。

那么如何控制混杂因素呢？——多因素分析：

多因素分析是相对于单因素分析而言，单因素分析仅关注一个因素在组间的差异或对结局事件的效应大小，而不考虑其他因素的影响。但实际上一种结局事件的发生和发展，常常受到多个因素的共同作用，因此仅采用单因素分析往往并不十分合理。多因素分析则是把多个变量之间的内在联系和相互影响考虑在内，同时分析多个因素对结局的影响。最常用的有3种：

在多因素回归分析中，不管是多重线性回归、logistic回归、还是Cox回归，通常的做法是，将我们在研究中关注的暴露/处理因素，以及可能的混杂因素一同放入到回归模型中进行拟合，如果模型显示暴露/处理因素对结局事件的效应值有统计学显著性，则可认为在“调整了”（Adjusted）其他混杂因素的影响后，该暴露/处理因素对于结局事件是一个“独立”（Independent）的影响因素。

这时就可以说该因素的独立危险因素了。

1.3.6 多因素怎么选择因素去做？

我们首先来看一篇2011年发表在The New England Journal of Medicine (影响因子：72.4)的文章：《A Prospective Natural-History Study of Coronary Atherosclerosis》：

Baseline variables that were considered clinically relevant or that showed a univariate relationship with outcome were entered into multivariate Cox proportional-hazards regression model. Variables for inclusion were carefully chosen, given the number of events available, to ensure parsimony of the final model.

重点有3个：

这些单因素应该和临床相关，如年龄、性别、吸烟等
从单因素的结果中筛选，但是考虑到单因素可能会被一些混杂因素所掩盖真正的结果，所以可以考虑将单因素p值放宽，在这篇新英格兰杂志中采用了单因素分析p<0.2。
根据样本数量考虑模型的稳健程度。具体如3.5 独立危险因素章节部分的图所示。

控制混杂因素的个数主要取決于发生结局事件的多少：

对于多重线性回旧模型，样本量应至少为10-15的自变量个数。也就是说，如果你打算做多重线性回归，想要在模型中纳入10个变量，那就要求样本量至少为100-150个。
对于logistic回归和Cox回归，结局事件则应至少为15-20倍的自变量个数。也就是说，如果你准备做 logistic或Cox回归，想要在模型中纳入10个变量，那么要求所需要的结局事件至少为150-200个。需要注意的是，这里指的是结局事件的数量，而不是总的样本量。

2. 用R语言进行Cox回归分析

强烈推荐——sthda官网R代码教程：http://www./english/wiki/cox-proportional-hazards-model

2.1 载入数据

library("survival")
library("survminer")
data("lung")
head(lung)

	inst	time	status	age	sex	ph.ecog	ph.karno	pat.karno	meal.cal	wt.loss
1	3	306	2	74	1	1	90	100	1175	NA
2	3	455	2	68	1	0	90	90	1225	15
3	3	1010	1	56	1	0	90	90	NA	15
4	5	210	2	57	1	1	90	60	1150	11
5	1	883	2	70	1	0	100	90	NA	0
6	12	1022	1	74	1	1	50	80	513	0

可以看到这里的因素不仅有分类变量（如性别），还有数值变量（如年龄，体重减轻）

2.2 单变量计算COX模型

直接使用函数即可：coxph(Surv(time, event) ~ variable, data, method)

# 单因素cox分析
if (T) {
  # 单个单因素分析
  head(lung)
  res.cox <- coxph(Surv(time, status) ~ ph.ecog, data = lung)
  summary(res.cox)

# 批量单因素分析
if (T) {
covariates <- c(“age”, “sex”, “ph.karno”, “ph.ecog”, “wt.loss”)
univ_formulas <- sapply(covariates,
function(x) as.formula(paste('Surv(time, status)~’, x)))

univ_models <- lapply( univ_formulas, function(x){coxph(x, data = lung)})
# Extract data
univ_results <- lapply(univ_models,
function(x){
x <- summary(x)

p.value<-signif(x$wald[“pvalue”], digits=2)
wald.test<-signif(x$wald[“test”], digits=2)

beta<-signif(x$coef[1], digits=2);#coeficient beta

HR <-signif(x$coef[2], digits=2);#exp(beta)
HR.confint.lower <- signif(x$conf.int[,”lower .95″], 2)
HR.confint.upper <- signif(x$conf.int[,”upper .95″],2)

HR <- paste0(HR, ” (“,
HR.confint.lower, “-“, HR.confint.upper, “)”)

res<-c(beta, HR, wald.test, p.value)
names(res)<-c(“beta”, “HR (95% CI for HR)”, “wald.test”,
“p.value”)

return(res)
#return(exp(cbind(coef(x),confint(x))))
})
res <- t(as.data.frame(univ_results, check.names = FALSE))
as.data.frame(res)
}
}

2.2.1 性别单因素结果分析

关于性别单因素分析的结果——res.cox <- coxph(Surv(time, status) ~ ph.ecog, data = lung)：

Call:
coxph(formula = Surv(time, status) ~ sex, data = lung)

n= 228, number of events= 165

coef exp(coef) se(coef)      z Pr(>|z|)
sex -0.5310    0.5880   0.1672 -3.176  0.00149 **
—
Signif. codes:  0 '***’ 0.001 '**’ 0.01 '*’ 0.05 '.’ 0.1 ' ’ 1

exp(coef) exp(-coef) lower .95 upper .95
sex 0.588 1.701 0.4237 0.816

Concordance= 0.579  (se = 0.021 )
Likelihood ratio test= 10.63  on 1 df,   p=0.001
Wald test            = 10.09  on 1 df,   p=0.001
Score (logrank) test = 10.33  on 1 df,   p=0.001

我们来理解下这个结果：

n= 228, number of events= 165 ——总人数是228个个体，有165个最后死了。
是否具有统计学意义：z给出的是Wald检验的统计量，z = coef/se(coef) 。而我们最关心的就在于，到底有没有意义呢？就看最后一列——Pr(>|z|)，因为这里的结果是0.00149，所以我们可以得出结论，性别差异具有统计学意义。
回归系数：coef这个值代表的是回归系数，具体是什么意思呢？就是说，这里我们用到的是以性别为分组信息，Male=1 Female=2 ，这里的coef=-0.5310，就是第二组相对于第一组来说的HR——女性比男性有更低的危险。
HR：`coef`的e次方结果，`exp(coef)`就是我们真正的HR结果了。这里，exp(coef)=0.588表示作为女性，危险性会下降58.8%，如果是男性的话，则会增加1.7倍危险性。
最后给出了一个模型的检验，有3种方法的结果，这3个意义相同。

2.2.2 批量单因素分析结果

放上批量单因素处理后的具体结果：

	beta	HR (95% CI for HR)	wald.test	p.value
age	0.019	1 (1-1)	4.1	0.042
sex	-0.53	0.59 (0.42-0.82)	10	0.0015
ph.karno	-0.016	0.98 (0.97-1)	7.9	0.005
ph.ecog	0.48	1.6 (1.3-2)	18	2.69e-05
wt.loss	0.0013	1 (0.99-1)	0.05	0.83

sex，age和ph.ecog变量具有较大的beta值，而ph.karno的beta则不显着。
age和ph.ecog具有正的β系数，而性别具有负β系数。因此，高龄和高phogog与较差的存活率有关，而女性（男1女 2）与较好的存活率有关。

2.3 多因素Cox分析

同样的函数调用，只不过后面需要多加几个变量值：

# 多因素cox分析
if (T) {
  res.cox <- coxph(Surv(time, status) ~ age + sex + ph.ecog, data =  lung)
  summary(res.cox)
}

2.3.1 多因素分析结果

Call:
coxph(formula = Surv(time, status) ~ age + sex + ph.ecog, data = lung)

n= 227, number of events= 164
(1 observation deleted due to missingness)

coef exp(coef)  se(coef)      z Pr(>|z|)
age      0.011067  1.011128  0.009267  1.194 0.232416
sex     -0.552612  0.575445  0.167739 -3.294 0.000986 ***
ph.ecog  0.463728  1.589991  0.113577  4.083 4.45e-05 ***
—
Signif. codes:  0 '***’ 0.001 '**’ 0.01 '*’ 0.05 '.’ 0.1 ' ’ 1

exp(coef) exp(-coef) lower .95 upper .95
age        1.0111     0.9890    0.9929    1.0297
sex        0.5754     1.7378    0.4142    0.7994
ph.ecog    1.5900     0.6289    1.2727    1.9864

Concordance= 0.637  (se = 0.025 )
Likelihood ratio test= 30.5  on 3 df,   p=1e-06
Wald test            = 29.93  on 3 df,   p=1e-06
Score (logrank) test = 30.5  on 3 df,   p=1e-0

结果的解读和上面的一样，但是我们可以发现这里的值和上面的不同了。这里我们可以探索下对于sex这个因素在单因素分析和多因素分析中的Cox结果：

	coef	HR
单因素	-0.5310	0.588
多因素	-0.552612	0.575445

可以看到在单因素和多因素Cox分析中，sex对于患者的影响都是一致的，女性会降低危险，但是具体大小会发生改变——因为多因素时考虑了其他因素的影响，但是coef值改变很小，说明性别对于生存有明显作用，且没有收到其他因素太大的影响。

2.4 可视化展示

得到模型后，我们想要对模型进行可视化展示，这时可以考虑使用森林图：

# 森林图可视化展示
if (T) {
  ggforest(res.cox,  #coxph得到的Cox回归结果
           data = lung,  #数据集
           main = 'Hazard ratio of data',  #标题
           cpositions = c(0.05, 0.15, 0.35),  #前三列距离
           fontsize = 1, #字体大小
           refLabel = 'reference', #相对变量的数值标签，也可改为1
           noDigits = 3 #保留HR值以及95%CI的小数位数
  )
}