TCGA系列学习笔记(6)多因素生存分析

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.mediecogroup.com/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.bioinfo-scrounger.com/archives/647/
  • Cox与KM生存分析及结果解读:https://www.omicsclass.com/article/1138
  • 强烈推荐——sthda官网R代码:http://www.sthda.com/english/wiki/cox-proportional-hazards-model
  • Cox回归生存分析 – 简书:https://www.jianshu.com/p/e80eb4168043
  • R语言实现及结果解读:https://www.omicsclass.com/article/1132
  • 【8文合集】全面了解单因素分析和多因素分析:https://www.mediecogroup.com/method_topic_article_detail/583/?ty=methods

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

https://www.mediecogroup.com/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 生存分析种类

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

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

备注

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

1.3 Cox回归

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

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

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

1.3.1 模型建立
  1. 不能直接以S(t,X)为因变量,以X为自变量进行线性回归的原因:
  • 生存分析研究中的数据包含删失数据
  • 时间变量t通常不满足正态分布和方差齐性
  1. Cox比例风险回归模型:以风险率函数h(t,X)为因变量,建立指数形式的回归方程(公式1):

 

h(t,X) = h0(t) exp( β1X1 + β2X2 + β3X3 + … + βmXm )
 

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

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

每个exp( βiXi ) 称为HR(hazard ratios):HR>1是危险因素;HR=1无影响;HR<1有利因素。

  1. Cox回归模型的优点
  • 未对h0(t)作任何假定,因此Cox回归模型在处理问题时具有较大的灵活性,属于非参数检验
  • 在许多情况下,即使h0(t)未知,仍可以估计出参数β。(意义:将时间和危险因素对风险函数的影响分开,仅仅研究危险因素对生存的影响,而不关心时间对生存的影响)
  1. 上述公式可以转化为广义线性模型,以便计算和解释(公式2):

 

ln[ h(t,X) / h0(t) ] = lnHR =  β1X1 + β2X2 + β3X3 + … + βmXm
 

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

1.3.2 基本假定
  • 比例风险假定

各危险因素的作用不随时间的变化而变化(公式3)

 

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

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

  • 对数线性假定

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

 

ln[ h(t,X) / h0(t) ] = lnHR =  β1X1 + β2X2 + β3X3 + … + βmXm
 

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

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

1.3.3 偏回归系数的意义

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

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

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

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

1.3.4 参数估计与假设检验
  • 其他自变量不变的情况下,变量 Xi 每增加一个单位,相对危险度RR的(1-α)置信区间为:

表示β的标准误。

  • 对于回归模型的假设检验通常采用似然比检验、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个:

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

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

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

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

强烈推荐——sthda官网R代码教程:http://www.sthda.com/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(&gt;|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
  • sexageph.ecog变量具有较大的beta值,而ph.karno的beta则不显着。
  • ageph.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的小数位数
  )
}

可以看到森林图左下角会给出出现结局事件的个数,COX生存模型的P值,AIC值和C-index值。

后记

后面应该就剩下最后的Cox模型前提假设的验证,以及对模型的表现进行打分,例如ROC,AUC了。后面再说吧~

生物信息学

TCGA系列学习笔记(4)差异分析结果注释

2020-10-7 22:22:17

生物信息学

TCGA系列学习笔记(7)建模及模型评价

2020-10-8 22:46:30

加入Q群
0 条回复 A文章作者 M管理员
    暂无讨论,说说你的看法吧
个人中心
今日签到
有新私信 私信列表
搜索