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

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

内容目录

前言1. 背景知识1.1 Cox前提假设的验证1.2 lasso和ridge回归1.3 C-index1.3.1 如何计算1.3.2 关于C-index的取值1.3.3 用R计算C-index1.4 Cox模型验证1.4.1 数据拆分1.4.2 ROC和AUC1.5 逐步回归1.6 nomogram列线图2. 代码教程2.1 Cox模型前提假设检验2.1.1 载入数据2.1.2 构建多因素模型2.1.3 模型前提假设检验——因素独立于时间2.1.4 可视化展示2.1.5 模型前提假设检验——对数风险比与危险因素线性相关2.2 Cox模型的建立流程2.2.1 载入数据2.2.2 批量单因素Cox分析2.2.3 多因素Cox分析2.2.4 逐步回归2.2.5 用建立好的模型进行预测2.2.6 时间依赖的ROC曲线2.2.6.1 survivalROC包2.2.6.2 timeROC包2.2.7 nomogram列线图后记

前言

这次的笔记其实在Cox模型的预测上卡了一段时间,不过好在咨询了曾老师后,得到了方向的指引,在此,由衷感谢曾老师的指点!

就用曾老板亲自编辑的感谢词来感谢吧:

Fat Yang thanks Dr.Jianming Zeng(University of Macau), and all the members of his bioinformatics team, biotrainee, for generously sharing their experience and codes !

1. 背景知识

1.1 Cox前提假设的验证

在之前的推文中讲到如果要使用cox模型分析,那么需要有2个关键的前提假设:

  • 各危险因素的作用不随时间的变化而变化

对于危险因素独立于时间的假设,从统计学上来说,我们可以基于scaled Schoenfeld residuals进行检验:

Schoenfeld residuals原则上独立于时间,如果发现Schoenfeld residuals与时间的关系不是随机的,p值显著,就说明模型违背了cox的假设

  • 对数风险比应与模型中的自变量应与呈线性关系

对于对数风险比与危险因素线性相关的假设,从统计学上来说,我们可以基于Martingale residuals进行检验:

通过绘制Martingale residuals 和连续性变量之间的关系图来观察。

Martingale residuals的取值范围在:-INF ~ +1

  • 接近1表示死的太快
  • 负值越大说明活得约久

注意:非线性相关的检验只对连续型变量进行检验

具体的代码见后面代码部分。

1.2 lasso和ridge回归

因为这里涉及到较多关于机器学习的概念,所以这节内容也是有点多的。

先来学习下基础的概念——lasso和ridge回归。

StatQuest – 正则化之岭回归_Ridge Regression

StatQuest – 正则化之 Lasso 回归

为什么会有lasso和ridge回归这2种回归呢?主要是为了解决机器学习过程中引起的过拟合现象。关于过拟合问题的解释可以去看:

回归问题-Lasso回归:https://blog.csdn.net/foneone/article/details/96576990

目前针对过拟合的问题,存在2种解决方案:

  • 方法一:减少特征数量(人工选择重要特征来保留,会丢弃部分信息)。
  • 方法二:正则化保留所有的特征变量,但是会减小特征变量的数量级。)。

lasso和ridge回归都是正则化的一种方式。

正则化的作用就是选择经验风险和模型复杂度同时较小的模型。

防止过拟合的原理:正则化项一般是模型复杂度的单调递增函数,而经验风险负责最小化误差,使模型偏差尽可能小经验风险越小,模型越复杂,正则化项的值越大。要使正则化项也很小,那么模型复杂程度受到限制,因此就能有效地防止过拟合。

因为目前我的数据不需要用到lasso回归,所以这部分代码先不学习了

1.3 C-index

https://www.bilibili.com/video/av68461486

C- index,英文名全称 concordance index,中文里有人翻译成一致性指数,最早是由范德堡大学( Vanderbilt University)生物统计教教授 Frank E Harrell Jr1996年提出,主要用于计算生存分析中的Cox模型预测值与真实之间的区分度( discrimination),和AUC其实是差不多的。

一般评价模型的好坏主要有两个方面

  • 一是模型的拟合优度(Goodness of Fit),常见的评价指标主要有R方、2logL、AIC、BlC
  • 另外一个是模型的预测精度,顾名思义就是模型的真实值与预测值之间差别大小,均方误差,相对误差等

在临床应用上更注重预测精度,建模的主要目的是用于预测,而C- index它就属于模型评价指标中的预测精度

1.3.1 如何计算
  1. 所有样本互相配对,共有N×(N-1)/2对,其中N为样本数,得到全部的对子数。

例如:A、B、C、D、E、F六个病人(N=6),可以配成6×(6-1)÷2=15对:

(A,B)(A,C)(A,D)(A,E)(A,F);

(B,C)(B,D)(B,E)(B,F);

(C, D)(C, E)(C, F);

(D,E)(D,F);

(E,F)

  1. 去除配对中两种情况
  • 两个病人都没有达到事件终点(比如死亡)
  • 其中的一个病人A的生存时间短于另一个病人B,然而病人A还没有到达事件终点(死亡)

因为这2种配对无法判断出谁先死的,此时剩下的配对数记为:M。

实际例子来看:

病人 生存时间 生存状态 描述 模型预测生存率
A 11 1 死亡 0.73
B 9 0 生存 0.69
C 14 1 死亡 0.76
D 3 1 死亡 0.70
E 4 0 生存 0.65
F 16 0 生存 0.79

那么这里需要剔除的有:

  • 情况1:需剔除(B,E)(B,F)(E,F)三对
  • 情况2:需剔除(A,B)(A,E)(B,C)(C,E)四对

拿其中一个来解释:(B,C)对中B的生存时间9个月短于C的生存时间14个月,B还未达到终点事件(依旧存活),这种情况下,因为B的随访时间短,和C的随访时间不一样长,不能判断B一定比C活的久,如果是相同随访时间下,B依旧活着,可以认为B的生存时间比C长。

剩下的对子数M=15-3-4=8对

  1. 计算剩下的配对中,预测结果和实际相一致的配对数记为K,即:两个病人如果真实生存时间较长的一位其预测生存时间长于另一位,或预测的生存概率高的一位的生存时间长于另位,则称之为预测结果与实际结果相符,称之为一致:

    例如:(AC)对子中,A活了11天,C活了14天,而模型结果中A是0.73,C是0.76,这样预测结果和真实结果保持了一致。相反,对于(BD)对子,B可以活≥9天,D只活了3天,而模型结果中B是0.69,D是0.70,也就是说模型认为D应该活的更久,这样预测结果和真实结果不相符。

    上述例子中剩下的8对(AC)(AD)(AF)(BD)(CD)(CF)(DE)(DF)中:

    (AC)(AD)(AF)(CD)(CF)(DF)六对预测结果与生存时间一致,即:K=6

  2. 计算C-index=K/M=6÷8=0.75
1.3.2 关于C-index的取值

从上述计算方法可以看出C- index:在0.5-1之间(任意配对随机情况下一致与不一致刚好是0.5的概率):

  • 0.5为完全不一致,说明该模型没有预测作用
  • 1为完全一致,说明该模型预测结果与实际完全一致

一般情况下C- index:

  • 在0.50-0.70为准确度较低
  • 在0.71-0.90之间为准确度中等
  • 高于090则为高准确度
1.3.3 用R计算C-index

其实从之前的推文中可以看到,在前面Cox模型的森林图中,已经显示了C-index的具体值,说明我们在用survival包的函数coxph计算模型的时候,已经得出了这个模型的C-index值,现在需要做的只是提取出来而已,因为比较简单,就直接把代码放在这里了:

library(survival) 
# 多因素Cox建模
res.cox <- coxph(Surv(time, status) ~ age + sex + ph.ecog, data =  lung)
sum.surv<- summary(res.cox)
# 结果提取
c_index <- sum.surv$concordance
c_index

1.4 Cox模型验证

参考:

一张图搞懂临床预测模型构建方法选择

其实目前的做法有2种:

  • 一个数据集做模型训练,一个数据集做模型验证。例如:TCGA训练模型,得到模型后用GEO做验证
  • 一个数据集经过数据拆分后,得到2个数据子集,其中1个做模型训练,另1个数据集做模型验证。
1.4.1 数据拆分

如果你的数据集样本够多,那么可以考虑把数据集进行拆分。在拆分的时候,我们可以简单考虑随机对半拆分,但是这会导致一些因素的分布不均匀,例如:500个样本进行随机拆分成2个250样本的数据集,数据集1中80%都是女性,数据集2中30%是女性,那么这就可能引起评价模型时的错误。

这时我们可以考虑使用caretR包,通过机器学习的方法帮我们得到最佳的切割分组,使得分组后样本因素分布的相对均匀。

函数说明:

Train_set_n <- createDataPartition(lung$status,p = 0.5,list = F)
# 常用选项
# y:样本的结局事件
# p:训练集占到的数量比例
# list:返回结果是否作为一个list对象
1.4.2 ROC和AUC

基于R软件实现联合诊断的ROC分析

在临床工作中,我们遇到的问题往往比较复杂,大多数疾病的诊断并非依靠单一诊断指标,我们往往需要检测不同指标后才能做出诊断。

比如,对于SLE的诊断,我们需要检测多个生化或免疫指标,并结合影像学检查与临床症状和体征做出诊断。实际临床实践中,联合诊断的情况更为常见。所谓联合诊断,是一种串联形式的诊断试验设计,即多个诊断指标阳性,患者患病的可能性更高。就和我们现在做的事情一样:为了找到一些可以指征预后的mRNA。

下面笔者提几个问题供大家思考:

(1). 什么情况下才需要联合不同的诊断指标?

(2). 如何确定联合多个诊断指标的诊断效能?

(3). 如何评价不同诊断指标联合方式的优劣?

针对第一个问题,如果单一指标的诊断效能不高,比如ROC曲线下面积AUC低于0.8,或者即便高于0.8,但临床有更高诊断效能的需求,都可以进行联合诊断

针对第二个问题,需要计算多个诊断指标联合的参数,一般是以参考标准的诊断结果为因变量,以待评价指标为自变量,构建Cox回归模型,并计算每个对象对应的预测概率,以预测概率进行ROC分析,计算曲线下面积AUC

针对第三个问题,比较不同联合诊断方式的ROC曲线下面积AUC即可。

1.5 逐步回归

https://blog.csdn.net/dingchenxixi/article/details/50543822

当我们做多因素分析的时候,问题来了,这时变量多了,该怎么选择合适的变量呢?

  • 选定变量(多元),因为有时候变量太多,有些对于模型是没有帮助甚至是倒忙
  • 多重共线性:有些变量可以由其他变量推出来
  • 检验模型是否合理

这里我们就可以使用到逐步回归——多元线性回归选择变量的方法:

  • 向前引入法:从一元回归开始,逐步增加变量,使指标达到最优
  • 向后剔除法:从全变量回归方程开始,逐步删去某个变量,使指标值达到最优为止
  • 逐步筛选法:综合上述两种方法

在变量选取之前,有几个判断指标先介绍一下

我们通过逐步回归,就是为了让模型在尝试各种变量组合的时候,使得AIC值最小,从而得到最佳的模型

具体算法细节可参考:https://wenku.baidu.com/view/0cd259ae69dc5022aaea0043.html

1.6 nomogram列线图

https://mp.weixin.qq.com/s/6vrWuOchNY_1qhs1pWBzrw

这个比较偏临床了,不过感觉画出来的图还是挺有趣的,所以也顺便整理下笔记吧!感觉这种图可以放到文章的最后一个图,还是挺不错的!后面就直接看代码的教程吧,在那里我再讲解下怎么去看这种图。

2. 代码教程

2.1 Cox模型前提假设检验

2.1.1 载入数据
# 载入数据及R包
if (T) {
  rm(list = ls())
  options(stringsAsFactors = F)
  library("survival")
  library("survminer")
  data("lung")
  head(lung)
}
2.1.2 构建多因素模型
# 构建多因素cox模型
if (T) {
  res.cox <- coxph(Surv(time, status) ~ age + sex + wt.loss, data =  lung)
  res.cox
}
2.1.3 模型前提假设检验——因素独立于时间
# 对模型进行检验
# 基于scaled Schoenfeld residuals进行检验
# Schoenfeld residuals原则上独立于时间,如果发现Schoenfeld residuals与时间的关系不是随机的,p值显著,就说明模型违背了cox的假设
if (T) {
  # 构建模型
  res.cox <- coxph(Surv(time, status) ~ age + sex + ph.ecog, data =  lung)
  summary(res.cox)

# 检验模型
test.ph <- cox.zph(res.cox)
test.ph
#          chisq df    p
# age     0.5077  1 0.48
# sex     2.5489  1 0.11
# wt.loss 0.0144  1 0.90
# GLOBAL  3.0051  3 0.39
}

可以看到对于每个因素来说p值都不显著,说明是符合cox的前提假设的,所以我们这个模型是可以用的。

2.1.4 可视化展示
# 可视化展示
if (T) {
  ggcoxzph(test.ph)
}

下图中,实线表示拟合的模型曲线,虚线表示mean ± 2sd 的置信区间:

如果发现有因素不独立于时间,那么可以:

  • 去除该因素
  • 根据因素分层进行检验
  • 加入一个时间*因素的interaction作用
2.1.5 模型前提假设检验——对数风险比与危险因素线性相关

因为只能对连续型变量进行分析,所以我们这里只检验age和wt.loss,考虑到wt.loss有1个数据缺失,所以我们将缺失的样本删除后再进行分析:

# 对模型进行检验
# 基于Martingale residuals进行检验
# 通过绘制Martingale residuals 和连续性变量之间的关系图来观察
if (T) {
  lung <- lung[!is.na(lung$wt.loss),]
  ggcoxfunctional(Surv(time, status) ~ age + wt.loss, data = lung)
}

可以看到年龄似乎还不错,是成正线性相关,而体重下降值似乎不那么明显,有些负线性相关的样子而已:

考虑到这里只是凭主观判断,没有p值选择,所以我们就容忍下。

2.2 Cox模型的建立流程

考虑到我自己的模型是Cox模型,所以,这里就放上关于Cox模型的分析方法。为了方便大家学习,所以我就从头到尾做一个完整的分析流程吧!

2.2.1 – 2.2.2都是之前已经写过的部分了,所以就很快带过了。

2.2.1 载入数据

# 载入数据
if (T) {
  rm(list = ls())
  options(stringsAsFactors = F)
  library(ROCR)
  library("survival")
  library("survminer")
  data("lung")
  head(lung)

# 考虑到NA会对分析有影响,所以这里先简单剔除
sapply(lung, function(x) sum(is.na(x)))
lung <- lung[,-c(1,9,10)]
lung <- lung[!(is.na(lung$ph.ecog) |
is.na(lung$ph.karno) |
is.na(lung$pat.karno)),]
}

2.2.2 批量单因素Cox分析

# 批量单因素分析
if (T) {
  covariates <- colnames(lung)[3:ncol(lung)]
  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)
}
##             beta HR (95% CI for HR) wald.test p.value
## age        0.018            1 (1-1)       3.6   0.057
## sex        -0.52    0.6 (0.43-0.83)       9.3  0.0023
## ph.ecog     0.47        1.6 (1.3-2)        17 4.3e-05
## ph.karno  -0.015      0.98 (0.97-1)       6.6    0.01
## pat.karno  -0.02   0.98 (0.97-0.99)        13 0.00039

2.2.3 多因素Cox分析

根据单因素Cox分析的结果,我们可以看到似乎单从P值来看都还不错。那么就先把所有的因素都纳入到多因素Cox模型中吧

# 建立多因素Cox模型
if (T) {
  fit <- coxph(Surv(time, status) ~ age + sex + ph.karno + ph.ecog, data =  lung)
  summary(fit)
}
## Call:
## coxph(formula = Surv(time, status) ~ age + sex + ph.karno + ph.ecog + 
##     pat.karno, data = lung)
## 
##   n= 223, number of events= 160 
## 
##                coef exp(coef)  se(coef)      z Pr(>|z|)   
## age        0.011383  1.011448  0.009510  1.197  0.23134   
## sex       -0.561464  0.570373  0.170689 -3.289  0.00100 **
## ph.karno   0.015853  1.015979  0.009853  1.609  0.10762   
## ph.ecog    0.565533  1.760386  0.186716  3.029  0.00245 **
## pat.karno -0.010111  0.989940  0.006881 -1.470  0.14169   
## ---
## Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
## 
##           exp(coef) exp(-coef) lower .95 upper .95
## age          1.0114     0.9887    0.9928     1.030
## sex          0.5704     1.7532    0.4082     0.797
## ph.karno     1.0160     0.9843    0.9965     1.036
## ph.ecog      1.7604     0.5681    1.2209     2.538
## pat.karno    0.9899     1.0102    0.9767     1.003
## 
## Concordance= 0.647  (se = 0.025 )
## Likelihood ratio test= 32.9  on 5 df,   p=4e-06
## Wald test            = 33  on 5 df,   p=4e-06
## Score (logrank) test = 33.79  on 5 df,   p=3e-06

可以看到似乎在多因素分析中,之前显著的因素如age 、ph.karno、pat.karno和age变得不显著了,可能是因为sex或者ph.ecog对其有混杂影响,而在纳入这些因素后,他们就不显著了

本来我们可以直接将ph.karnoage直接剔除模型,但在实际分析过程中,可能经过这些操作后还会有很多基因存在,这时我们就需要进行逐步回归了,下面进行逐步回归检验操作。

2.2.4 逐步回归

# 逐步回归
if (T) {
  step.multi_COX= step(fit,direction = "both")
  # 综合两个结果,考虑保留sex和ph.ecog
  fit <- coxph(Surv(time, status) ~ sex + ph.ecog, data =  lung)
}
## Start:  AIC=1422.96
## Surv(time, status) ~ age + sex + ph.karno + ph.ecog + pat.karno
## 
##             Df    AIC
## - age        1 1422.4
## <none>         1423.0
## - pat.karno  1 1423.1
## - ph.karno   1 1423.7
## - ph.ecog    1 1430.3
## - sex        1 1432.4
## 
## Step:  AIC=1422.42
## Surv(time, status) ~ sex + ph.karno + ph.ecog + pat.karno
## 
##             Df    AIC
## <none>         1422.4
## - pat.karno  1 1422.6
## - ph.karno   1 1422.6
## + age        1 1423.0
## - ph.ecog    1 1429.7
## - sex        1 1431.7

可以看到逐步回归的过程是先在age + sex + ph.karno + ph.ecog + pat.karno5个变量的模型中进行分析:

  • 原始模型的AIC:1423.0
  • 去除age的AIC:1422.4
  • 去除pat.karno的AIC:1423.1
  • 去除ph.karno的AIC:1423.7
  • 去除ph.ecog的AIC:1430.3
  • 去除sex的AIC:1432.4

可以看到在去除age后的AIC最小,所以接下来去除age,得到的模型再进行分析逐步回归分析:

  • 原始模型的AIC:1422.4
  • 去除pat.karno的AIC:1422.6
  • 去除ph.karno的AIC:1422.6
  • 加上age的AIC:1423.0
  • 去除ph.ecog的AIC:1429.7
  • 去除sex的AIC:1431.7

这时发现原始模型的AIC最小,所以终止分析,得到最终模型是Surv(time, status) ~ sex + ph.karno + ph.ecog + pat.karno

但是我们检查模型会发现ph.karno、ph.karno、和pat.karno是不显著的,所以手动去除。

2.2.5 用建立好的模型进行预测

这里就需要了解一个函数——predict

具体参考:

  • predict函数中type参数的选择:https://bbs.pinggu.org/thread-6394043-1-1.html

上面的链接引入了一个问题:Cox模型由于baseline的危险值不知道,所以没有办法像其他模型那样,根据一个新的表达矩阵预测得到生存状态,所以如何预测Cox模型的好坏呢?下面的解决方案

  • Predictions for Cox model:http://personal.psu.edu/drh20/R/html/survival/html/predict.coxph.html
  • 如何验证Cox模型的AUC:https://www.biostars.org/p/109215/
  • predict.coxph结果的理解:https://stats.stackexchange.com/questions/44896/how-to-interpret-the-output-of-predict-coxph
  • 还可以参考这个文章:https://www.sciencedirect.com/science/article/pii/S0014483518306031?via%3Dihub
  • 或者直接在R中 ?predict.coxph 也可以看到帮助文档。

最重要的就是要明白:

  1. 由于Cox模型本身的建模特点,所以我们无法直接预测患者的生存或者死亡概率,因为我们不知道baseline的危险值
  2. 但是通过predict函数,我们可以通过type参数对Cox模型预测得到不同的数据,例如:
  • “lp”——当所有危险因素都用均值表示时的logHR的值.
  • “risk”——得到的是危险分数HR,也就是exp(lp).
  • “expected”——可能出现的死亡人数.
# 用构建的模型进行预测
# 自己预测自己
if (T) {
  # RiskScore<-predict(fit,type = "risk",newdata = lung[,c('sex','ph.ecog')])
  RiskScore<-predict(fit,type = "risk",newdata = lung[,c('sex','ph.karno','ph.ecog','pat.karno')])
  risk_group<-ifelse(RiskScore>=median(RiskScore),'high','low')

# 生存曲线
ggsurvplot(survfit(Surv(time, status) ~ risk_group, data=lung),
pval = TRUE, #show p-value of log-rank test,显示log-rank分析得到的P值
conf.int = TRUE, #添加置信区间
conf.int.style = “step”,  # customize style of confidence intervals,改变置信区间的样子
risk.table = “abs_pct”,  # absolute number and percentage at risk,这里以n(%)的形式展示risk table
risk.table.y.text.col = T,# colour risk table text annotations.
risk.table.y.text = FALSE,# show bars instead of names in text annotations in legend of risk table.不显示注释名字
xlab = “Time in days”,   # customize X axis label.自定义x的标签为time in years
surv.median.line = “hv”, #添加中位生存时间的线
ncensor.plot = FALSE, #我这里不显示删失的图,TRUE就显示
legend.labs = c(“high risk”, “low risk”),    # 对legend的标签重新命名
palette = c(“#E7B800”, “#2E9FDF”), # 自定义颜色
ggtheme = theme_light()#绘图主题
)
}

根据预测出来的HR进行分组后,对样本进行KM生存曲线的绘制,结果如下:

似乎还可以!不过也是情理之中,毕竟用的是同一套数据做出来的模型。不过没有出现死亡的全在1组,生存的全在一组,说明没有过拟合

那么接下来就去测评下这个模型的ROC曲线了。

2.2.6 时间依赖的ROC曲线

所谓时间依赖的ROC曲线,和普通ROC曲线不同的在于它可以预测在指定之间内,模型的ROC曲线。简单来说,就是可以预测一个模型1年生存率、3年生存率等不同时间点下的准确情况。

这里我们用2个R包去做。

2.2.6.1 survivalROC包
# 方法1——survivalROC
if (T) {
  library(survivalROC)
  survival_ROC<-survivalROC(Stime=lung$time, #生存时间,Event time or censoring time for subjects
                            status=lung$status-1, #生存状态,dead or alive
                            marker=RiskScore, #风险得分,Predictor or marker value    
                            predict.time=250, #预测250天的生存时间
                            method="KM" #使用KM法进行拟合,默认的方法是method="NNE" 
  )
  survival_ROC_AUC <- round(survival_ROC$AUC,3)
  survival_ROC_AUC

#画图
#x轴为False positive rate,y轴为True positive rate
plot(survival_ROC$FP,survival_ROC$TP,type=”l”,xlim=c(0,1),ylim=c(0,1),
xlab=”False positive rate”,
ylab=”True positive rate”,
main=paste0(“ROC Curve”, ” (“, “AUC = “,survival_ROC_AUC,” )”),  #标题
cex.main=1.5,#放大标题
cex.lab=1.3,#坐标轴标签(名称)的缩放倍数
cex.axis=1.3, font=1.2, #放大轴上的数字
lwd=1.5, #指定线条宽度
col=”red”  #红色
)
abline(a=0,b=1) #添加一条y=x的线
}

这里预测了250天的ROC曲线:

结果似乎不太好?难道模型没有建好?再用另外一个试试。

2.2.6.2 timeROC包
# 方法2——timeROC
if (T) {
  library(timeROC)
  time_ROC<-timeROC(T=lung$time, #生存时间(dead和alive的生存时间).
                    delta=lung$status-1, #生存结局,Censored的样本必须用0表示
                    marker=RiskScore, #预测的变量,这里是风险评分,在没有特殊情况下,值越大,越容易发生危险事件
                    cause=1, #阳性结局的赋值(必须是1或2),也就是dead的赋值,这里dead是1表示的
                    weighting = "marginal", #权重计算方法,这是默认方法,选择KM计算删失分布,weighting="aalen" [选用COX],weighting="aalen" [选用Aalen]
                    times = c(250,500,750), #计算250,500,750的ROC曲线
                    ROC=TRUE,
                    iid=TRUE #计算AUC
  )
  #绘制250的ROC
  plot(time_ROC,time=250,col="red",title=FALSE,lwd=2) #将生成一条两倍于 默认宽度的线条 
  #在此基础上添加5年的ROC
  plot(time_ROC,time=500,col="blue",add=TRUE,title=FALSE,lwd=2) 
  #add 10年的ROC
  plot(time_ROC,time=750,col="black",add=TRUE,title=FALSE,lwd=2)
  #添加图例
  legend("bottomright", #图例设置在右下角
         c(paste0("AUC at 250 days = ", round(time_ROC$AUC[1],3)),
           paste0("AUC at 500 days = ", round(time_ROC$AUC[2],3)), 
           paste0("AUC at 750 days = ", round(time_ROC$AUC[3],3))),
         col=c("red","blue","black"),lwd=2,bty="n" #或者bty+“o"
  )
}

这里求了3个时间的ROC曲线:

但是250天的结果和上面的一致。都不太行。

难道真的是模型的问题,于是我使用逐步回归得到的模型,再次尝试,结果如下:

方法1:

方法2:

看来还是用逐步回归得到的模型效果更好!

关于代码,还是一样的,就贴在最后以便以后回顾吧:

# 载入数据
if (T) {
  rm(list = ls())
  options(stringsAsFactors = F)
  library(ROCR)
  library("survival")
  library("survminer")
  data("lung")
  head(lung)

# 考虑到NA会对分析有影响,所以这里先简单剔除
sapply(lung, function(x) sum(is.na(x)))
lung <- lung[,-c(1,9,10)]
lung <- lung[!(is.na(lung$ph.ecog) |
is.na(lung$ph.karno) |
is.na(lung$pat.karno)),]
}

# 批量单因素分析
if (T) {
covariates <- colnames(lung)[3:ncol(lung)]
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)
}

# 建立多因素Cox模型
if (T) {
fit <- coxph(Surv(time, status) ~ age + sex + ph.karno + ph.ecog + pat.karno, data =  lung)
summary(fit)
}

# 逐步回归
if (T) {
step.multi_COX= step(fit,direction = “both”)
# 综合两个结果,考虑保留sex和ph.ecog
# fit <- coxph(Surv(time, status) ~ sex + ph.ecog, data =  lung)
fit <- step.multi_COX
}

# 用构建的模型进行预测
# 自己预测自己
if (T) {
# RiskScore<-predict(fit,type = “risk”,newdata = lung[,c(‘sex’,’ph.ecog’)])
RiskScore<-predict(fit,type = “risk”,newdata = lung[,c(‘sex’,’ph.karno’,’ph.ecog’,’pat.karno’)])
risk_group<-ifelse(RiskScore>=median(RiskScore),’high’,’low’)

# 生存曲线
ggsurvplot(survfit(Surv(time, status) ~ risk_group, data=lung),
pval = TRUE, #show p-value of log-rank test,显示log-rank分析得到的P值
conf.int = TRUE, #添加置信区间
conf.int.style = “step”,  # customize style of confidence intervals,改变置信区间的样子
risk.table = “abs_pct”,  # absolute number and percentage at risk,这里以n(%)的形式展示risk table
risk.table.y.text.col = T,# colour risk table text annotations.
risk.table.y.text = FALSE,# show bars instead of names in text annotations in legend of risk table.不显示注释名字
xlab = “Time in days”,   # customize X axis label.自定义x的标签为time in years
surv.median.line = “hv”, #添加中位生存时间的线
ncensor.plot = FALSE, #我这里不显示删失的图,TRUE就显示
legend.labs = c(“high risk”, “low risk”),    # 对legend的标签重新命名
palette = c(“#E7B800”, “#2E9FDF”), # 自定义颜色
ggtheme = theme_light()#绘图主题
)
}

# 时间依赖的ROC曲线
if (T) {
# 方法1——survivalROC
if (T) {
library(survivalROC)
survival_ROC<-survivalROC(Stime=lung$time, #生存时间,Event time or censoring time for subjects
status=lung$status-1, #生存状态,dead or alive
marker=RiskScore, #风险得分,Predictor or marker value
predict.time=250, #预测250的生存时间
method=”KM” #使用KM法进行拟合,默认的方法是method=”NNE”
)
survival_ROC_AUC <- round(survival_ROC$AUC,3)
survival_ROC_AUC

#画图
#x轴为False positive rate,y轴为True positive rate
plot(survival_ROC$FP,survival_ROC$TP,type=”l”,xlim=c(0,1),ylim=c(0,1),
xlab=”False positive rate”,
ylab=”True positive rate”,
main=paste0(“ROC Curve”, ” (“, “AUC = “,survival_ROC_AUC,” )”),  #标题
cex.main=1.5,#放大标题
cex.lab=1.3,#坐标轴标签(名称)的缩放倍数
cex.axis=1.3, font=1.2, #放大轴上的数字
lwd=1.5, #指定线条宽度
col=”red”  #红色
)
abline(a=0,b=1) #添加一条y=x的线
}

# 方法2——timeROC
if (T) {
library(timeROC)
time_ROC<-timeROC(T=lung$time, #生存时间(dead和alive的生存时间).
delta=lung$status-1, #生存结局,Censored的样本必须用0表示
marker=RiskScore, #预测的变量,这里是风险评分,在没有特殊情况下,值越大,越容易发生危险事件
cause=1, #阳性结局的赋值(必须是1或2),也就是dead的赋值,这里dead是1表示的
weighting = “marginal”, #权重计算方法,这是默认方法,选择KM计算删失分布,weighting=”aalen” [选用COX],weighting=”aalen” [选用Aalen]
times = c(250,500,750), #计算250,500,750的ROC曲线
ROC=TRUE,
iid=TRUE #计算AUC
)
#绘制250的ROC
plot(time_ROC,time=250,col=”red”,title=FALSE,lwd=2) #将生成一条两倍于 默认宽度的线条
#在此基础上添加500的ROC
plot(time_ROC,time=500,col=”blue”,add=TRUE,title=FALSE,lwd=2)
#add 750的ROC
plot(time_ROC,time=750,col=”black”,add=TRUE,title=FALSE,lwd=2)
#添加图例
legend(“bottomright”, #图例设置在右下角
c(paste0(“AUC at 250 days = “, round(time_ROC$AUC[1],3)),
paste0(“AUC at 500 days = “, round(time_ROC$AUC[2],3)),
paste0(“AUC at 750 days = “, round(time_ROC$AUC[3],3))),
col=c(“red”,”blue”,”black”),lwd=2,bty=”n” #或者bty+“o”
)
}
}

2.2.7 nomogram列线图

代码其实比较简单,主要介绍下结果怎么解读吧

放上代码和结果:

# Nomogram图
if (T) {
  library(rms)
  # 打包数据
  if (T) {
    dd=datadist(lung)
    options(datadist="dd") 
  }

# 构建符合要求的模型
if (T) {
f1 <- psm(Surv(time,status) ~ age+ph.karno+ph.ecog+pat.karno, data =  lung, dist=’lognormal’)
med <- Quantile(f1) # 计算中位生存时间
surv <- Survival(f1) # 构建生存概率函数
}

# 绘制COX回归中位生存时间的Nomogram图
if (T) {
nom <- nomogram(f1, fun=function(x) med(lp=x),
funlabel=”Median Survival Time”)
plot(nom)
}

## 绘制COX回归生存概率的Nomogram图
## 注意lung数据的time是以”天“为单位
if (T) {
nom <- nomogram(f1, fun=list(function(x) surv(365, x),
function(x) surv(730, x)),
funlabel=c(“1-year Survival Probability”,
“2-year Survival Probability”))
plot(nom, xfrac=.6)
}
}

这里画了2种图:

  • 中位生存时间的Nomogram图

  • 生存概率的Nomogram图

其实可以把它画的更好看,只是鉴于时间紧张,今天就先这样了,不过放上一个链接,有兴趣的可以去学习下如何画出更好看的Nomogram图

基于R绘制竞争风险模型列线图

后记

今天内容有点多了,耐心慢慢看吧~

在此还是要感谢下曾老师,因为当我筛出基因后,因为建模是Cox模型,而Cox模型的ROC没法直接用predict预测生死,一度卡住了,经过曾老师指导后选择用预测得到的危险得分分组绘制KM生存曲线图,用预测得到的危险得分去绘制ROC曲线,一下就解决了我的问题!

就用曾老板亲自编辑的感谢词来感谢吧:

Fat Yang thanks Dr.Jianming Zeng(University of Macau), and all the members of his bioinformatics team, biotrainee, for generously sharing their experience and codes !

生物信息学

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

2020-10-7 22:22:50

基础实验

玩转qPCR系列:数据分析第二关

2020-8-14 22:27:05

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