简介面向数据分析人员与临床科研工作者的风险预测建模资料包围绕ROC曲线、PR曲线、NRI重分类等指标提供完整的模型构建与评估方案。压缩包共40个文件大小584KB以28个MATLAB脚本为主涵盖AUC比较、NRI计算、风险分层绘图等核心功能另含4个xlsx数据表、2个xls预实验数据及说明文档便于复现与二次开发。资料还包含无类别限制NRI、多类别NRI等进阶实现可辅助处理二分类与多分类预测场景适用于医疗预后、金融信贷等领域的风险概率估算。已有494人学习下载适合需要快速搭建预测模型并进行性能验证的研究者可结合自带数据与脚本直接开展实测。1. 523例预实验样本里的风险预测模型不是画一条ROC那么简单拿到“风险预测模型1225.rar”时我以为只是一份随手画一条ROC曲线的Matlab脚本。解压后看到里面躺着523例预实验数据、一堆NRI.m、AUC_compare_correlated.m和Risk_Assessment_Plot.m才意识到这套代码的价值不在画图而是把“模型好不好”这件事从单一AUC扩展到了重分类改善、多分类NRI、相关AUC比较的完整验证链。对于做临床预测模型或金融违约评分的从业者这套RAR包里的函数至少能帮你省掉三周的调试时间。下面我会从文件结构、ROC/PR计算到NRI与AUC比较一步步还原这套风险预测模型的可复现路径。2. 解析RAR包结构与数据入口main.m、zhen.m与Excel数据读取风险预测模型的交付方式通常是RAR包里面散落着函数脚本和Excel数据。放在这个包里的523例预实验(1).xls和风险20181225.xlsx是模型的输入源main.m和zhen.m则负责把数据装配成后续指标脚本能吃的格式。理解这个结构比急着跑main.m更重要。标题末尾的will7jv更像是发布者留下的批次标识代码本身不依赖这个字段但项目路径里最好别留这串字符免得Matlab导入时把文件当成函数名解析。2.1 解压RAR包与文件分类在Linux环境下解压$ unrar x 风险预测模型1225.rar如果机器上没装unrar可以用7-zip的命令行版本$ 7za x 风险预测模型1225.rarunrar x里的x参数表示保留压缩包内完整路径直接在当前目录展开7za同理。Windows下建议用WinRAR解压到纯英文路径例如D:\risk_model\免得Matlab因为中文路径读不出Excel。解压后还需要清理Mac归档元数据$ rm -rf __MACOSX $ find . -name ._* -delete这两行命令适用于任何来源的RAR包。__MACOSX和._*.m是macOS Finder生成的隐藏文件在Windows下虽然不显示但Matlab的addpath(genpath(.))会把空白的._NRI.m当成同名函数加载导致真正调用的NRI.m被覆盖。先清理再进Matlab能避免一类特别隐蔽的报错。解压后文件大致可以归成下面几类类别文件名预期作用主脚本main.m, zhen.m数据装配与流程串联指标计算Roc.m, CIAUC.m, NRI.m, multi_category_NRI.m, multi_category_NRI_ci.mROC坐标、AUC置信区间、重分类改善指标模型比较AUC_compare_correlated.m, Category_Free_NRI.m, Category_Free_NRI_ci.m相关AUC比较与连续NRI图表Risk_Assessment_Plot.m, ROC test.m风险分层图、ROC图数据523例预实验(1).xls, 风险20181225.xlsx, roc.xlsx, test.xlsx, origin.xlsx训练/验证数据、中间结果工具函数exciseRows.m, FindNandD.m, choice.m, LR_Pz_choice.m行清洗、事件数查找、逐步回归变量选择系统残留__MACOSX, ._*.m, ~$ADE ME RAP...rtfmacOS归档元数据和Office锁文件~$开头的是Office正在编辑该文件时留下的锁文件可以忽略。真正要关心的只有表格里前五行以及Excel数据列名与脚本之间的对应关系。2.2 main.m的数据装配Excel读取与列对齐大多数临床风险预测模型的main.m里前二十行都是同一件事读Excel、找标签列、找风险评分列。由于原包里的Excel是从统计软件导出的列名和数据类型不一定干净直接用readtable的前置探测参数会比较稳% 风险预测模型入口读取523例预实验数据 opts detectImportOptions(523例预实验(1).xls); opts.VariableNamingRule preserve; data readtable(523例预实验(1).xls, opts); % 用列名定位标签和风险评分不依赖硬编码列号 event data{:, event}; risk data{:, risk}; fprintf(有效样本%d事件数%d事件发生比例%.2f%%\n, ... height(data), sum(event), 100*mean(event));detectImportOptions会自动识别文件分隔符和字段类型VariableNamingRule设为preserve后Excel里的中文列名不会被Matlab翻译成Var1这样你才能放心地用data{:, event}。如果同一份数据里有事件列但没有专门的risk列就得先确认哪几列是特征再对特征做逻辑回归预测概率——这个包里的LR_Pz_choice.m干的就是这件事按P值筛选变量生成最终风险评分。老版本Matlab如果读不了.xlsx可以退回到xlsread[~, ~, raw] xlsread(风险20181225.xlsx);但xlsread会丢精度数字读回来变成数值数组字符串列变成cell还需要手工映射列名。新版本能用readtable就别走回头路。2.3 zhen.m的“中间人”角色关于zhen.m最合理的作用是组装多个模型的风险概率。假设你训练了三个模型逻辑回归、随机森林、XGBoost它们各自的预测概率散布在不同变量里后续Roc.m和NRI.m都要求输入等长的列向量。zhen.m就是把它们拼成一个矩阵% 把多模型预测概率拼成矩阵方便批量比较 riskMatrix [risk_logistic, risk_rf, risk_xgb]; labels event; % 逐列计算ROC for i 1:size(riskMatrix, 2) [x, y, ~, auc] perfcurve(labels, riskMatrix(:, i), 1); fprintf(模型%d AUC%.3f\n, i, auc); end这里perfcurve是Matlab自带的二分类评价函数第一个参数是真实标签第二个是预测概率第三个1说明正类标签是数值1。安全起见运行前先看一眼每个模型的预测概率是否都在0到1之间如果出现NaN或负数会导致画图时曲线断掉。zhen.m还有一种可能是把数据按分层抽样切分训练集和验证集特别是523例样本量不大时五折交叉验证比固定切分更可靠。这个包的exciseRows.m应该是配合行删除用的比如删除事件缺失的行。运行时可以先看FindNandD.m里有没有对事件数量的判断它通常在清洗步骤之后根据正负类样本量决定是否继续后续计算。3. ROC与PR曲线计算原理在Matlab中复现Roc.m与Risk_Assessment_Plot.mRAR包里的Roc.m是核心画图脚本Risk_Assessment_Plot.m则是风险分层图。画图本身不难难的是理解为什么同一份数据会产生ROC和PR两种看起来打架的曲线。3.1 从混淆矩阵到TPR、FPR的定义二分类模型输出一个预测概率后我们会选一个截断阈值。高于阈值为“预测阳性”低于阈值为“预测阴性”。由此得到四个格子TP实际阳性预测阳性、FP实际阴性预测阳性、TN实际阴性预测阴性、FN实际阳性预测阴性。真阳性率 TPR TP / (TP FN)反映所有真实阳性里有多少被识别出来假阳性率 FPR FP / (FP TN)反映所有真实阴性里有多少被误判成阳性。ROC曲线就是在变阈值时把每个(FPR, TPR)点连起来。PR曲线则是横轴召回率Recall等同TPR纵轴精确率Precision TP / (TP FP)。它的关键区别是ROC把阴性样本的分母算进FPR里当阴性样本非常多时FPR对误判不敏感PR曲线把精确率作为纵轴假阳性一多精确率立刻跳水。3.2 用perfcurve同时得到ROC与PR曲线perfcurve不仅画ROC也能画PR。关键在于xCrit和yCrit两个参数% 利用perfcurve实现ROC与PR的对称绘制 labels event; scores risk; % ROC曲线XFPRYTPR [rocX, rocY, ~, aucRoc] perfcurve(labels, scores, 1); % PR曲线X召回率Y精确率 [prX, prY] perfcurve(labels, scores, 1, ... xCrit, reca, yCrit, prec); figure(Position, [100 100 900 380]); subplot(1, 2, 1); plot(rocX, rocY, b-, LineWidth, 2); hold on; plot([0 1], [0 1], k--); xlabel(False Positive Rate); ylabel(True Positive Rate); title(sprintf(ROC (AUC%.3f), aucRoc)); xlim([0 1]); ylim([0 1]); grid on; subplot(1, 2, 2); plot(prX, prY, r-, LineWidth, 2); xlabel(Recall); ylabel(Precision); title(Precision-Recall Curve); grid on;对于ROC模式perfcurve的第一个输出是FPR第二个输出是TPR对于PR模式第一个输出是召回率第二个输出是精确率。AUC计算结果只在ROC模式下有意义PR模式下的“曲线下面积”受插值方式影响不建议直接当作最终指标。如果需要给ROC加置信区间可以追加NBoot参数% 用1000次Bootstrap给ROC加置信区间 [~, ~, ~, aucRoc, ~] perfcurve(labels, scores, 1, ... NBoot, 1000, BootType, normal);BootType可选normal、per、bca。风险预测模型里样本量仅523例我一般选normal它计算速度最快且结果稳定bca更准确但Bootstrap次数不足时容易失败。参数说明如下perfcurve参数作用常见取值labels真实标签0/1必须只有两个取值否则用 unique 检查scores模型预测概率连续值允许单调变换posclass指定哪个类别为正类1 或 positivexCrit横轴准则tpr默认ROCrecaPRyCrit纵轴准则fpr默认ROCprecPRNBoot自助法置信区间默认0需要置信区间时设1000BootType置信区间类型normal / per / bca3.3 Risk_Assessment_Plot.m背后的分层风险图Risk_Assessment_Plot.m在实际项目中承载的是“预测风险 vs 观测风险”的校准图。最常见的实现是把样本按预测概率从低到高分五组然后看每组里实际事件比例是否递增。这比单条ROC更能说明模型能不能直接用于分层决策% 按预测概率分5层计算每层实际事件率 edges linspace(0, 1, 6); riskGroup discretize(risk, edges); observed zeros(1, 5); for g 1:5 sel (riskGroup g); if sum(sel) 0 observed(g) mean(event(sel)); end end bar(1:5, observed, FaceColor, [0.3 0.6 0.8]); xlabel(Predicted Risk Group (1Low, 5High)); ylabel(Observed Event Rate);这里linspace(0, 1, 6)把0到1等切成5段discretize把每个风险概率归入对应分组。要注意的是如果样本集中在0.1附近固定阈值分箱会让中间组出现空档画出来的图会有零柱。更稳妥的做法是用prctile(risk, 0:20:100)按分位数分组保证每组样本量大致相等。3.4 ROC输出异常时先查这三处临床数据画ROC最常出现的三个问题是第一标签列不是严格0/1而是字符串“yes/no”perfcurve会报错第二预测概率列里有缺失值直接导致曲线断点第三正类选错方向画出来的曲线在对角线下方。我在复现这套代码时习惯在画图前加一行unique(labels)验证标签取值再用isnan(risk)删除缺失样本。这些校验逻辑都放在代码前面编译时不会被注意到但能避免一大批低级错误。4. NRI重分类改善与AUC比较multi_category_NRI.m和AUC_compare_correlated.m的使用边界RAR包里价值最高的函数我认为不是Roc.m而是NRI.m、multi_category_NRI.m和AUC_compare_correlated.m。在风险预测模型里常见场景是老模型加了新标记物后AUC只提升了0.01但把低危患者正确重分类到了高危组或者把高危患者正确降到了低危组。这种改善在ROC上几乎看不到却直接影响临床决策。4.1 从分类表到NRI公式净重分类改善指数统计的是相比旧模型新模型是否把更多事件组患者正确升到更高风险层以及更多非事件组患者正确降到更低风险层。令“向上移动”表示新模型比旧模型把样本分到更高风险层“向下移动”相反。NRI的计算公式对于事件组实际发生事件的人NRI_event P(up | event) - P(down | event)对于非事件组NRI_nonevent P(down | nonevent) - P(up | nonevent)总NRI NRI_event NRI_nonevent。正的NRI意味着新模型整体重分类正确率更高。要注意的是这个定义依赖风险分层的阈值。比如把 10% 视为低危10%~20% 视为中危20% 视为高危分层的区间不同NRI数值也会不同。Category_Free_NRI.m是连续NRI不需要人为设定分层用样本在全部分位点上的移动汇总multi_category_NRI.m则支持三个及以上风险类别。NRI类型适用场景是否需要指定阈值分类NRI临床上有明确风险分层界值是连续NRI探索性分析、还没有公认分层标准否多分类NRI风险等级超过两档是每级都要界值4.2 在原代码中调用NRI函数的最小姿势由于原包里的NRI.m我们没法直接看到内部签名但按照这类工具包的惯例函数输入大概率是“事件标签、旧模型预测概率、新模型预测概率”。调用时先用 try-catch 包裹出错了就读代码头注释% 用try-catch包裹NRI调用避免函数签名不符时中断 event data.event; risk_old data.old_risk; risk_new data.new_risk; try [nri_point, p_value] NRI(event, risk_old, risk_new); fprintf(分类NRI %.3f, p %.4f\n, nri_point, p_value); catch ME if strcmp(ME.identifier, MATLAB:undefinedFunction) disp(NRI.m不存在或文件不在当前路径); else rethrow(ME); end end这里nri_point是点估计p_value是正态近似的显著性检验。NRI的分布不像AUC那样容易用大样本公式所以原包里的multi_category_NRI_ci.m应当是用Bootstrap给NRI算置信区间。我们可以用同样的思路自己补一个1000次自助法% 自助法计算NRI的95%置信区间 nBoot 1000; n numel(event); bootNRI zeros(nBoot, 1); rng(2025); for b 1:nBoot idx randsample(n, n, true); try bootNRI(b) NRI(event(idx), risk_old(idx), risk_new(idx)); catch bootNRI(b) NaN; end end nriCI quantile(bootNRI(~isnan(bootNRI)), [0.025, 0.975]); fprintf(NRI 95%% CI [%.3f, %.3f]\n, nriCI(1), nriCI(2));randsample(n, n, true)表示从1到n中有放回地抽取n个样本每次采样样本量与原始数据集相同。Bootstrap的置信区间如果严重跨零说明两模型的改善不稳定即使点估计为正论文里也不能下“有显著改善”的结论。4.3 相关AUC比较与NRI的配合AUC_compare_correlated.m解决的是另一个问题两个模型都用同一份数据训练它们预测概率之间有相关性直接对AUC做t检验是错的需要用DeLong检验或Bootstrap配对比较。下面这段代码演示了用perfcurve每次取出AUC再做配对Bootstrap% 配对Bootstrap比较两个模型的AUC差值 nBoot 1000; aucDiff zeros(nBoot, 1); rng(42); for b 1:nBoot idx randsample(n, n, true); [~, ~, ~, aucOld] perfcurve(event(idx), risk_old(idx), 1); [~, ~, ~, aucNew] perfcurve(event(idx), risk_new(idx), 1); aucDiff(b) aucNew - aucOld; end pAUC mean(aucDiff 0) * 2; % 双侧检验近似如果aucDiff的分布绝大多数都大于0说明新模型AUC确实更高。这个结果和NRI放一起看能区分“整体区分度提升”和“重分类方向正确”两种不同的改善避免只看单一指标被误导。4.4 多分类NRI的使用边界multi_category_NRI.m把二分类的NRI推广到多分类。对于高危、中危、低危三个类别事件组的移动就不只是向上或向下而是从真实类别向预测类别的偏移总和。这套口径计算复杂度高样本量不足时很容易产生奇异值。523例样本做二分类NRI问题不大强行套多分类NRI必须保证每个风险层都有足够事件数——一般经验是最小层事件数不少于30。否则我宁愿退回连续NRI也决不用分类NRI硬撑。5. 把模型迁移到自己的数据从替换Excel到校验ROC输出的完整套路5.1 替换数据的三步核对把这套RAR包迁移到自己的数据最重要的不是改代码而是确认三个前提第一标签必须是0/1数值列不能有字符串第二模型预测概率列必须与标签行一一对应删过行之后要重新排序对齐第三Excel里不能有合并单元格readtable遇到合并单元格会把左侧值填为空。核对完成后把523例预实验(1).xls换成你自己的数据文件名main.m里的列名改成你的列名。原包里的origin.xlsx和test.xlsx很可能就是作者留下的“原始数据”和“验证数据”两种切分建议沿用这种命名方便回溯。5.2 校验输出的简易hook我习惯在Roc.m调用结束后用下面的代码做一次独立校验% 递推法验证ROC曲线坐标的单调性 [rocX, rocY] perfcurve(event, risk, 1); assert(numel(rocX) 10, ROC点太少阈值变化不足); assert(all(diff(rocX) 0), FPR非单调递增检查标签方向); assert(all(diff(rocY) 0), TPR非单调递增检查排序); fprintf(ROC坐标校验通过点数%d\n, numel(rocX));如果断言报错优先检查标签方向是否反了再检查预测概率是否重复值过多。正常模型的ROC坐标点来自每个不重复概率阈值点数过少意味着评分被硬四舍五入成了整数值。5.3 验证数据文件与脚本输出的交叉对照原RAR包里的roc.xlsx我拿到手的第一反应是拿它和Roc.m输出做减法。先读文件% 读取作者上次保存的ROC结果并对比 refROC readtable(roc.xlsx); diffX abs(refROC.x - rocX); if max(diffX) 0.01 warning(与ROC参考结果偏差超过0.01检查数据版本); end这种对照不需要跑完整模型却能快速发现上一手数据是否被清洗过。偏差大的时候回到main.m检查exciseRows.m是否删除了违规行。5.4 解压时容易掉进的三个坑这个RAR包解压后第一坑是__MACOSX文件夹混在路径里Matlab的addpath(genpath(.))会把里面的无用m文件也加进来而有些._NRI.m是空白文件导致调用NRI时加载到空函数。第二坑是~$开头的Office锁文件它不影响Matlab运行但如果你用dir(*.m)批量递归添加路径文件名会混入非法字符。第三坑是Excel文件被另一个Excel进程占用readtable会报权限错误关掉Excel再跑一次即可。样本量敏感时最后做一步把roc.xlsx输出与自己的ROC结果做差偏差超过0.01就说明某处阈值或标签方向不一样。把这个核对函数和5.2的hook合并成一个check_roc_output.m放进公共脚本库以后任何二分类模型在换数据后都能直接复用。本文还有配套的精品资源点击获取
