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 如何计算

所有样本互相配对，共有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）

去除配对中两种情况：

两个病人都没有达到事件终点（比如死亡）
其中的一个病人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对

计算剩下的配对中，预测结果和实际相一致的配对数记为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
计算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$concordancec_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.karno和age直接剔除模型，但在实际分析过程中，可能经过这些操作后还会有很多基因存在，这时我们就需要进行逐步回归了，下面进行逐步回归检验操作。

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 也可以看到帮助文档。

最重要的就是要明白：

由于Cox模型本身的建模特点，所以我们无法直接预测患者的生存或者死亡概率，因为我们不知道baseline的危险值
但是通过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——survivalROCif (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——timeROCif (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)

# 批量单因素分析
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’)

# 时间依赖的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

# 方法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列线图

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

怎么解读：https://www.mediecogroup.com/method_topic_article_detail/371/?ty=methods
怎么画：基于R绘制竞争风险模型列线图

放上代码和结果：

# 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 !

本站仅提供存储服务，所有内容均由用户发布，如发现有害或侵权内容，请点击举报。