预测模型 | 4. Nomogram(列线图)的构建、美化及其风险分数的计算原理
目 录
一、Nomogram构建前的准备
二、Nomogram代码详解
三、Nomogram美化
四、Nomogram风险分数的计算原理
五、代码、练习数据
1. 理清思路
2. 变量转换
R语言 | 18. 连续变量取最佳cutoff值: 基于KM曲线; R语言 | 19. 连续变量取最佳cutoff值: 基于ROC曲线; R语言 | 20. 连续变量取最佳cutoff值: 基于限制性立方样条 (RCS)
3. 变量筛选
载入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. 构建多因素回归模型
我们已经通过三种筛选变量的方法比较和确认了 年龄(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)
由于我们在前面已经将变量变为正确的英文格式,因此构建好的列线图的英文名字或缩写以及数轴上的亚变量的名字都是标准的。无需再用修图软件修改。
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强制标记的每个刻度线都绘制标签,我也没研究明白
)
#运行以展示多因素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)
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。转化为风险分数:
nomo2#运行后将得出所有变量的分数
其计算结果与我们手动计算的完全一致,此结果也向我们展示了每个刻度的分数为0.02429555 ,它正好是最大回归系数的百分之一。
end