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

目  录

一、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   npredict(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 

手动计算风险分数

1、模型的回归系数β(coef)

回归系数展示了某变量对模型的贡献,其绝对值越大贡献越大,正负值表示对结局的影响为正向还是反向。
同时,风险比HR值也是由回归系数算出来的,公式为HR=exp(coef),coef=回归系数。

2. Nomogram的风险分数是基于回归系数计算的

所有变量的回归系数都是以自己的某个亚变量为参考计算的(设参考变量coef=0)。

例如:rt:yes=-1.2162,是相对于No(其coef设为0)来说,Yes会使复发风险下降。但是这个参考不能适用于淋巴结转移n。Nomogram是一个整体,必须要有一个所有变量都能参考的值。

本代码里的列线图参考是最大回归系数。即拥有最大的回归系数的变量n=≥4,coef=2.4296在列线图里的分数为为100分。其他所有变量以它为标准进行转化。

3. 分类变量风险分数计算

以g变量作为例子:组织学分级II级相对于I级(coef=0)回归系数=1.3497。则转换公式为:

组织学分级为II级在Nomogram上的分数为55.55≈56,与图中一致

4. 连续变量的风险分数计算

首先,连续变量的回归系数解释为每变动1个单位其风险的变化。因此,计算出1个单位的分数就能算出所有分数。

以年龄为例:coef=-0.0315,说明每减小增加1岁风险减小0.0315。转化为风险分数:

结果:年龄变化1个单位,即1岁的风险分数=1.2963分
若80岁时患者复发风险为0分,则20岁时为1.2963*(80-20)=77.79≈78分
nomo2#运行后将得出所有变量的分数

其计算结果与我们手动计算的完全一致,此结果也向我们展示了每个刻度的分数为0.02429555 ,它正好是最大回归系数的百分之一

这就是Nomogram的本质,即使你不会R代码,也可以通过计算在word或PPT中画出Nomogram。

end

(0)

相关推荐