预测模型 | 4. Nomogram（列线图）的构建、美化及其风险分数的计算原理

阿非ycfg 2021-04-20

展开全文

目录

一、Nomogram构建前的准备

二、Nomogram代码详解

三、Nomogram美化

四、Nomogram风险分数的计算原理

五、代码、练习数据

一、构建Nomogram前的准备

1. 理清思路

特别注意，事先阅读既往文章，了解需要收集的变量和结局等注意事项，点击：

2. 变量转换

1-连续变量是否要转为分类变量，截断值怎么选；

2-有序的分类变量是否应合并等（如淋巴结转移按N1-3分还是按是否有转移分类），点击：

3. 变量筛选

最佳的变量组合应根据临床意义而定，而非统计学指标，然而，大多数文章是参考统计学指标，点击：

二、构建Nomogram的代码

载入R包和数据

#1.加载R包
library(rms)

#2.载入数据，status=0为复发
aa<- read.csv('nomo示例.CSV')
#3. 查看数据类型，本例结局为status，0为复发
str(aa)

本例数据涉及毕业论文不予公开，但是有类似的练习数据供大家下载使用。
str( )函数查看数据的性质，除time和age其他应均为分类变量。
同时，变量及其亚变量应变为正规的英文名字以方便投稿使用。

1、转换变量性质及其名字

#数值变量转为因子，
#且不管原始数据里用abc还是012表示，均转为其真正的名字。
#列线图里的文字将以修改后的文字显示
a$n<-factor(aa$n, 
             levels=c('a','b','c'),  
             labels=c('0','1-3','≥4')) 
aa$hr<- factor(aa$hr, 
               levels=c('a','b'),  
               labels=c('Negative ','Positive')) 
aa$lvi<-factor(aa$lvi, 
               levels=c('0','76'),  
               labels=c('Absent','Present'))
aa$g<-factor(aa$g, 
             levels=c(1,2,3), 
             labels=c('I','II','III'))
aa$rt<-factor(aa$rt, 
              levels=c('ano','shi'), 
              labels=c('No','Yes'))

factor()函数是将数值转为因子的意思。
例：aa数据的g变量（组织学分级），其亚变量为123，默认为连续变量，factor（aa$g）将其转为分类变量，且函数将123，赋值为I,II,II。这样列线图会显示I,II,III

#再次查看数据性质
str(aa)

结果：除time和age保留为连续变量（int），其他变量均转为分类变量（factor），且亚变量也变为英文

#变量名字改为其本来的英文意思
label(aa$age)<-'Age at diagnosis'
label(aa$n) <-'LNM'
label(aa$hr)<-'HR status'
label(aa$g)<-'Histologic grade'
label(aa$lvi)<-'Lymphovascular invasion'
label(aa$rt)<-'Radiation Therapy'

2. 构建多因素回归模型

在推文R语言 | 预测模型变量筛选: 代码篇中
我们已经通过三种筛选变量的方法比较和确认了 年龄(age)，淋巴结转移数（n），激素受体（hr），组织学分级（g），脉管癌栓（lvi），放疗（rt）为最佳变量组合。因此，我们这里不用再选择，可直接使用这5个变量来建模。

#加载数据，
#rms包下面两行代码需要执行，不然易报错。
nomo<-datadist(aa)
options(datadist='nomo')

#构建多因素回归模型
nomo1 <- cph(Surv(time,status==0)~age+n+er+lvi+g,
             x=T,y=T,
             data=aa,
             surv=T,
             time.inc = 12*5);nomo1

#运行结果
#                      Model Tests       Discrimination    
# Indexes        
# Obs        388    LR chi2    152.14    R2       0.372    
# Events      72    d.f.            8    Dxy      0.716    
# Center -0.2539    Pr(> chi2) 0.0000    g        1.546    
# Score chi2 207.28    gr       4.691    
# Pr(> chi2) 0.0000                      
# 
#             Coef    S.E.    Wald Z  Pr(>|Z|)
# age         -0.0315 0.0108  -2.90   0.0037  
# n=1-3        1.3234 0.5028   2.63   0.0085  
# n=≥4         2.4296 0.4730   5.14  <0.0001 
# hr=Positive -1.0217 0.2770  -3.69   0.0002  
# lvi=Present  1.3742 0.2932   4.69  <0.0001 
# g=II         1.3497 0.3386   3.99  <0.0001 
# g=III        1.3940 0.3822   3.65   0.0003  
# rt=Yes      -1.2162 0.2673  -4.55  <0.0001

代码解释：
示例数据结局：time=随访时间（月），status=复发否，==0为复发。【若你的数据1为复发或死亡，则改为==1】。
time.inc=12*5，现在没啥用，这里在之后的校准图要对应起来。
其他地方与多因素cox回归一致。
注意rms的Cox回归函数名为cph( )
结果解释：
结果汇报了每个变量的p值，回归系数β等
R²：评价模型拟合情况
Dxy：计算C指数：C-index=Dxy/2+0.5
回归系数β（coef）：列线图中的分数以此为基础计算，很重要！！

#或者直接计算C指数
Cindex <- rcorrcens(Surv(as.numeric(aa$time),aa$status==0)~predict(nomo1))
Cindex

                   C    Dxy  aDxy    SD     Z P   n
predict(nomo1) 0.142 -0.716 0.716 0.045 15.78 0 388

C-index =1-c； =aDxy/2+0.5
se=sd/2；
标准误=标准差/2；
95%CI=c-index +/- 1.96*se ；
C-index约等于ROC曲线的AUC值。
若使用coxph()函数，它会直接显示c指数和标准误se

3.设置预测时间

#设置预测时间
#里预测3个：1,3,5年
surv <- Survival(nomo1)#使用刚才构建的模型预测
surv1 <- function(x)surv(12*1,lp=x)
surv2 <- function(x)surv(12*3,lp=x)
surv3 <- function(x)surv(12*5,lp=x)

4. 构建列线图

nomo2<-nomogram(nomo1,
                #1.感兴趣的预测时间
                fun=list(surv1,surv2,surv3),
                #2. 预测时间段的名字
                funlabel=c('1-year LR probability',
                           '3-year LR probability',
                           '5-year LR probability'),
                #3.是否显示回归系数轴，一般为F
                lp =T, 
                #4.#分数为百分制
                maxscale=100,
                #5.设置预测复发率的范围【根据自己数据修改】
                fun.at=c('0.99','0.95','0.9','0.8',
                         '0.7','0.6','0.5','0.4',
                         '0.3','0.2','0.1')
                );plot(nomo2)

由于我们在前面已经将变量变为正确的英文格式，因此构建好的列线图的英文名字或缩写以及数轴上的亚变量的名字都是标准的。无需再用修图软件修改。

三、Nomogram美化

plot(nomo2, 
     
     #1.变量与图形的占比
     xfrac=.35, 
     #2.变量字体加粗
     cex.var=1, 
     #3.数轴：字体的大小
     cex.axis=0.8,
     #4.数轴：刻度的长度
     tcl=-0.5,
     #5.数轴：文字与刻度的距离
     lmgp=0.3, 
     #6.数轴：刻度下的文字，1=连续显示，2=隔一个显示一个
     label.every=1, 
     #7.1个页面有几个数轴(这个可以压缩行间距)
     naxes=13,
     #8.垂直线的颜色.
     col.grid=gray(c(0.8, 0.95)),
     #9.线性预测轴名字
     lplabel='Linear Predictorlp', 
     #10变量分数名字
     points.label='Points', 
     #11总分名字
     total.points.label='Total Points',
     force.label=T#没啥用。TRUE强制标记的每个刻度线都绘制标签，我也没研究明白
     )

四、Nomogram变量的赋值原理

构建的多因素Cox回归时，我们得出了每个变量的回归系数β（即coef）；在构建的Nomogram中，我们也可以得到每个变量的具体分数。

预测模型本质上就是多因素Cox回归的可视化而已，那Nomogram的风险分数是如何计算的呢？

就像森林图，它展示的是Cox回归的HR及95%CI，而Nomogram是将回归系数标准化，然后以数轴上以风险分数显示。

#运行以展示多因素COX回归结果的回归系数
nomo1
#运行以展示Nomogram每个变量的分数
nomo2

             Coef    S.E.   Wald Z Pr(>|Z|)
 age         -0.0315 0.0108 -2.90  0.0037  
 n=1-3        1.3234 0.5028  2.63  0.0085  
 n=≥4         2.4296 0.4730  5.14  <0.0001 
 hr=Positive -1.0217 0.2770 -3.69  0.0002  
 lvi=Present  1.3742 0.2932  4.69  <0.0001 
 g=II         1.3497 0.3386  3.99  <0.0001 
 g=III        1.3940 0.3822  3.65  0.0003  
 rt=Yes      -1.2162 0.2673 -4.55  <0.0001