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

风险预测模型评价之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就是当time2000，且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, fam