风险预测模型评价之NRI的R语言计算

风险预测模型评价之NRI的R语言计算R语言代码的盛（暴）宴（击），除了NRI的运算，还有蛮多预处理的操作，在哪都能用得着。大家做好战斗准备。R里有2个包专事计算NRI，分别为nricens和PredictABEL。从最后结果来说，nricens计算出来的是绝对NRI，PredictABEL则为相对NRI。但我们已经知道计算原理了呀，而且它们都能生成新旧模型分类的对照表，所以其实只用其中一个包就都可以计算了。不过它们还是有一些小小差异，我们就以logistic回归模型为例，分别看一下这两个包，供大家参考选择。Cox模型参数较多也较复杂，但我相信你看完这篇的讲解就能看懂帮助文档中的cox案例，算是留个小作业给你吧~拟合模型先用一份示例数据做个模型。这是survival包里带有的一份梅奥诊所的数据，记录了418位患者的临床指标，观察这些因素与原发性胆汁性肝硬化（PBC）的关系。其中前312位参加了RCT，其他的只参加了观察队列。我们用前312份样本，做个预测2000天时间点上死亡与否的模型。先载入这份数据看一下。library(survival)

###logistic回归

egData<-pbc[1:312,]做一个logistic回归，我们需要一个结局事件作为因变量，它必需是个分类变量；其次需要若干自变量，它们可以是分类也可以是连续。这个表中的结局是status，0=截尾（删失），1=接受移植，2=死亡。研究目的“死亡与否”是个二分类变量，所以要做些变换。再看time一栏，有的不够2000天，这些样本要么是没到2000天就死亡了，要么是删失了。我们要删掉2000天内删失的数据。egData=egData[egData$time>2000|

(egData$time<2000&egData$status==2),]“[]”表示筛选条件，|表示“或”，&为“和”。所以条件句就是egData中的time一列大于2000的保留，或小于2000但同时状态为死亡的也保留。最后一个“,”别忘了，其在条件句的前面表示对列进行选择，在其后表示选择行。选好后做一个event向量，把status的三个状态变成死亡=1,未死亡=0。event=ifelse(egData$time<2000&egData$status==2,1,0)ifelse(test,yes,no)大法好啊，前面一个test是逻辑判断句，其值为真时返回yes的值，为假时返回no的值。所以本句中test就是当time<2000，且status为2（死亡）时，记为1，否则为0。然后把event合并入原来的表格。egData=cbind(egData,event)cbind()是以列合并，另有rbind()以行合并。这样event就成了最后一列，为结局事件。然后选择模型的自变量（predictors）。太多了，选取其中几个示例。就以年龄、胆红素、白蛋白为旧模型（standard），三者加上一个凝血酶原时间为新模型（new）。一般做logistic回归是用glm(因变量~自变量1+自变量2+……+自变量n,

family=binomial('logit'),

data=数据表)，但如果自变量较多的话，前面那个运算式就会很长很长，万一这些自变量还是基因名或编号，就很想死了。所以顺便讲一个简化的办法，即把那一串先写成formula。fml.std<-as.formula(paste('event~',

paste(colnames(egData)[c(5,11,13)],

collapse='+')))这里有好几层函数，paste()会把括号中的元素粘贴起来，collapse是其中的间隔。colnames()是获取表格的列名，[]中的数值向量为所选择的列序号。这样如果是一个超大表格，你选中第10~70列还可以写成“10:70”。好了，同样写出新模型的formula：fml.new<-as.formula(paste('event~',

paste(colnames(egData)[c(5,11,13,19)],

collapse='+')))

可以查看一下，新模型的formula写成这个效果：然后像上边说的那样用glm()拟合两个模型。mstd=glm(fml.std,family=binomial('logit'),

data=egData,x=TRUE)

mnew=glm(fml.new,family=binomial('logit'),

data=egData,x=TRUE)

这样一长串运算式用刚才命名好的fml.std和fml.new代替就好了。x=TRUE是将来用nricens包计算时要求用到的，表示输出结果中是否包含所用到的数据表，平时可以不写。模型就这样做完了~先不急着计算NRI，先看看它的总体情况。summary(mstd)运行这句就得到该模型的描述特征。残差、相关系数、各个自变量的统计显著性等，注意倒数第二行的AIC，就是上一期提到的赤池信息准则，表示模型校准度，很少有人汇报呢。可以用同样的方法看看新模型。这里就-不展开了，进入下一环节。NRI的计算•先看nricens包的方法。library(nricens)NRI<-nribin(mdl.std=mstd,mdl.new=mnew,

updown='category',cut=c(0.3,0.6),

niter=10000,alpha=0.05)

填上新旧两个模型。Cut是判断风险高低的临界值，现在我们写了2个，也就是0~29%为低风险，30%~59%为中风险，60%~100%为高风险。现实中可以查阅相关文献进行设置，预测风险达到多少需要怎样干预之类的。Updown为定义一个样本的风险是否变动的方式，category是指分类值，即我就熟悉的低、中、高风险，另有一种diff，为连续值。选diff时，cut就设1个值，比如0.02，即认为当预测的风险在新旧模型中相差2%时，即被认为是重新分类了。这种方法用的比较少。后面几个参数就比较有意思了，niter为重复取样的次数，即boostrap方法，不做的话将其设为0就好了；做的话建议至少1000次，这也是默认值，但我（读书少）见过的研究都10000次。然后将统计显著性alpha设为0.05。这样就可以看到输出的结果：如果不做bootstrap，就是这个结果。有重新分类情况的详表，最后是NRI和各种变动的概率。第一个NRI如前所述，是绝对NRI，大家可以根据之前的知识和上边的详表自己计算验证一下，此时可手动计算出相对NRI。其他指标随便看看。如果做了bootstrap，就会多出一个表：因为做了10000次重复取样，相当于有10000个NRI，于是就有了标准误和置信区间，刚才我们设alpha=0.05，所以后面的Lower和Upper就是95%置信区间的下界和上界。同时，做不做bootstrap都会得到一张图，表示各数据点在新旧模型中的分布。默认的Case和Control标签我觉得不太严谨，Case代表结局事件中编号为“1”的组，也就是发生了结局（死亡），Control为“0”，未发生。其实是positive和negative比较贴切吧。反正就这个意思。这张图也和重新分类表的意思差不多，看看就好。•再看PredictABEL包的做法library(PredictABEL)

pstd<-mstd$fitted.values

pnew<-mnew$fitted.values

先把两个模型中的预测风险值提出来，也就是模型中的fitted.value。这个包只能从预测风险计算，刚才的nricens包可以用模型，也可以用预测风险（把mdl.std和mdl.new参数换成p.std和p.new）。reclassification(data=egData,cOutcome=21,

predrisk1=pstd,predrisk2=pnew,

cutoff=c(0,0.30,0.60,1))cOutcome是结局事件的列序号，刚才我们不是把event放到最后了么，即第21列，填上。两个预测风险值也相应填上。这里的cutoff跟刚才的不一样，还要填上前面的0和后面的1，成为完整的0~100%的区间。然后得到一个重新分类表：跟上边nricens做的差不多了。不过这个包没有bootstrap的选项。接着看下面的结果，这里的分类NRI是咱们上回说的相加NRI，同样可以根据上一期的知识手动计算一下。记得咱们并没有设置bootstrap吧？可这里也有个95%置信区间，只是内部调用了一个更为简陋的只能计算连续NRI的improveProb()函数做的，而且连续变化的临界值也不太透明，遂不管了。最后还有个IDI是指，发生和未发生结局事件样本的平均预测