简介:面向医学科研与金融风控场景的MATLAB风险预测模型工具包,涵盖ROC与PR曲线、NRI重分类、AUC比较等核心评估方法,适用于临床预测模型构建、疾病复发或违约风险评分等二分类任务,并可处理样本不平衡下的模型性能评估。压缩包共40个文件,其中28个.m文件为主要算法脚本,4个xlsx/xls为示例数据与结果表,3个rtf为说明文档,另有txt说明与脚本备份,整个压缩包仅584KB,便于快速下载与部署。目前已有498人学习下载使用。使用者可直接运行Roc.m、main.m等脚本,借助内置数据复现523例预实验的分析流程,并通过Risk_Assessment_Plot模块输出可视化风险图表;多类别NRI与AUC相关性比较工具则在新旧指标对比、模型优化中提供现成实现,适合需要快速搭建预测模型评估体系、开展风险量化研究的研究者与数据分析师。
1. 从523例样本到ROC曲线:这套风险预测模型代码包到底解决了什么
先说一个反直觉的结论:这套风险预测模型代码包里最值得研究的不是Roc.m绘图脚本,而是数据清洗和变量筛选那部分——因为最终ROC预测结果能不能站得住,关键在523例样本进入模型之前发生了什么。这是一份用MATLAB写成的风险预测模型实现,覆盖了从Excel原始表到逻辑回归建模、ROC曲线绘制、AUC置信区间与模型对比、NRI净重分类改善指数计算以及风险评估图输出的完整链路。压缩包里的origin.xlsx和test.xlsx是已经切分好的训练与验证数据。适合正在做二分类预测、手头有临床或金融表格数据、希望直接照着改并快速拿到AUC与ROC结论的研究者。
2. 数据准备阶段:523例样本进入模型前的三道清洗关口
2.1 先把压缩包里的文件分清楚:哪些是核心脚本、哪些是数据、哪些可以直接删
打开"风险预测模型1225.rar"后,我习惯先做一件事:把文件分成四组。第一组是数据文件,包括523例预实验.xls、origin.xlsx、test.xlsx、roc.xlsx;第二组是数据清洗与变量选择的脚本,包括exciseRows.m、zhen.m、choice.m;第三组是建模与评估脚本,包括LR_Pz_choice.m、Roc.m、ROC test.m、CIAUC.m、AUC_compare_correlated.m;第四组是NRI相关脚本,包括NRI.m、Category_Free_NRI.m、Category_Free_NRI_ci.m、multi_category_NRI.m、multi_category_NRI_ci.m、FindNandD.m。压缩包名字里的1225是打包日期,will7jv应该是上传者标识,和代码逻辑没有关系,不影响运行。
把分组说得更细一点。zhen.m从命名看是"真实"的拼音,对应的是真实标签处理;FindNandD.m是统计事件组与非事件组人数的前置函数,为NRI计算做准备。另外压缩包根目录下有__MACOSX文件夹,以及一堆以._开头的文件,这些是macOS系统自动生成的元数据,跟分析没有关系。还有"~$"开头的rtf文件,是Office打开文档时留下的临时文件。这些东西在Windows上解压后不会自动消失,建议解压后先手动清掉,否则MATLAB扫描当前目录时会把隐藏文件也扫进去,偶尔出现莫名其妙的路径报错。我把这份文件清单整理成表格,方便你对照使用:
| 分组 | 文件 | 作用 |
|---|---|---|
| 数据 | 523例预实验.xls、origin.xlsx、test.xlsx、roc.xlsx | 训练/验证原始表与结果导出 |
| 清洗 | exciseRows.m、zhen.m、choice.m | 行剔除、标签对齐、变量选择 |
| 建模 | LR_Pz_choice.m | 逻辑回归P值逐步筛选 |
| ROC | Roc.m、ROC test.m、CIAUC.m、AUC_compare_correlated.m | ROC绘制、置信区间、模型比较 |
| NRI | NRI.m、Category_Free_NRI.m、multi_category_NRI.m等 | 重分类指标计算 |
| 可视化 | Risk_Assessment_Plot.m | 风险分层图输出 |
把这张表放在手边,后面每一步要动哪个文件,直接对号入座就行。
2.2 行剔除与变量选择:exciseRows.m和choice.m的实际用法
exciseRows.m解决的是最基础的脏活:按行索引剔除不符合纳入标准的样本。
% exciseRows.m —— 按行索引剔除指定样本,返回清洗后的数据矩阵 % 输入: data为m×n数值矩阵,rowsToExclude为待剔除的行索引向量 function dataClean = exciseRows(data, rowsToExclude) if isempty(rowsToExclude) dataClean = data; return; end keepIdx = true(size(data, 1), 1); % 初始化全量保留掩码 keepIdx(rowsToExclude) = false; % 将待剔除行位置置为false dataClean = data(keepIdx, :); % 用逻辑索引一次性取出保留行 end逻辑说明:这里用逻辑掩码而不是data(rowsToExclude, :) = [],是因为后者在循环里多次删除时会导致行号整体前移,后一次删除的行号已经不是原始表里的行号了。用逻辑索引直接生成一个"保留哪些行"的真值数组,清洗完还能对照keepIdx检查到底剔除了哪几条样本,出问题时有后悔药可吃。
参数说明:rowsToExclude可以是行号向量,也可以是逻辑向量。我的经验是先从原始Excel里把不符合条件的人在"Excel行号"记下来,但要注意读入MATLAB后行号会少一行(表头占了第一行),所以传参之前先执行rowsToExclude = rowsToExclude - 1这种偏移,再核对一遍再跑,能少踩一个很隐蔽的坑。
choice.m则是按变量名筛选要进入建模的列:
% choice.m —— 按变量名返回需要建模的列子矩阵与列索引 % 输入: X为n×p数据矩阵,varNames为p个变量名的元胞数组, % selectedNames为要选出的变量名列表 function [Xsub, colIdx] = choice(X, varNames, selectedNames) [~, colIdx] = ismember(selectedNames, varNames); if any(colIdx == 0) error('变量名不存在: %s', strjoin(selectedNames(colIdx == 0), ', ')); end Xsub = X(:, colIdx); end逻辑说明:使用ismember做"按名引用"而不是硬编码列序号,是因为当从523例预实验.xls换成自己的数据时,Excel列顺序大概率会变,但列名通常是稳定的。这个脚本的价值在于:一旦你写错变量名,它会直接报错而不是默默返回一个空矩阵,这种显式失败比静默错误友好得多。
参数说明:varNames是从Excel表头读入的变量名元胞数组;selectedNames是你在分析计划里确定要纳入的预测因子。返回值里colIdx是这些变量在原矩阵中的实际列号,画相关性矩阵或做共线性检查时会经常用到。
2.3 Excel数据读入阶段的三个高发坑:文本型数字、缺失编码和列对齐
数据读入阶段最常见的问题不是统计上的,而是Excel本身带来的。第一坑是文本型数字:一列"血压值"如果部分单元格是文本格式,xlsread读进来后这一列会被切分成两列,或者直接变成NaN。第二坑是缺失值编码不统一:有人用空白、有人用0、有人输入"NA",如果缺失被编码成0,后续模型会把0当成真实测量值,整个预测概率分布都会偏移。第三坑是列对齐:在Excel里做过vlookup或者手动粘贴的人都知道,某一列与其他列的行序经常对不齐,这种事用代码很难自动发现。
我处理这些表时有一个固定顺序:先用readtable整表读入,对每一列统计NaN比例和唯一值个数;然后把所有字符串列统一转成double,并记录转换失败的位置;最后对关键变量做散点图矩阵,用肉眼确认行序没有错位。这套动作很基础,但已经帮我避免过三次"AUC很高,后来发现标签错位"的事故。记住一个原则:所有高级统计指标的前提,是表干净。
3. 逻辑回归变量筛选:LR_Pz_choice.m的P值门槛与训练/验证数据切分
3.1 为什么风险预测模型优先选择逻辑回归而不是机器学习
风险预测场景的结局大多是二分类:疾病发生或未发生、违约或不违约。逻辑回归输出的是事件概率,天然可以映射成风险评分;每个变量的系数能翻译成OR值(优势比),配合95%置信区间,这是临床论文和业务评审里最熟悉的语言。换成随机森林或XGBoost,预测精度可能更高,但可解释性弱很多,参数调整路径也长,在只有523例的小样本环境下还容易过拟合。所以我的习惯始终是:先用逻辑回归建一个基线模型,拿到AUC和校准度指标后,再决定要不要上更复杂的模型。这套代码包里的LR_Pz_choice.m,做的就是逻辑回归里最关键的一步——变量筛选。
3.2 LR_Pz_choice.m的工作逻辑:P值门槛怎么设、循环怎么停
LR_Pz_choice.m从名字看就是"Logistic Regression P value Choice",核心思路是逐步回归,用P值决定哪些变量进来、哪些出去。在压缩包里,它和单独跑一个glmfit的区别在于:它是自动化的、可复现的变量筛选过程,而不是手动一个一个试变量。我用伪代码还原它的主循环,方便你理解每一步在干什么:
% LR_Pz_choice.m —— 基于P值的逐步逻辑回归变量筛选(主循环示意) % 输入: X为n×p预测因子矩阵,y为n×1二分类结局(0/1) % pEnter: 变量进入模型的P值上限,默认0.05 % pRemove: 变量被剔除的P值下限,默认0.10 function [selected, coef, stats] = LR_Pz_choice(X, varNames, y, pEnter, pRemove) if nargin < 4, pEnter = 0.05; end if nargin < 5, pRemove = 0.10; end selected = []; % 当前入选变量集合 candidateIdx = 1:size(X, 2); % 尚未进入模型的候选变量 while true % 第一步:对每个候选变量单独做单因素logistic回归,记录P值 % 第二步:找出最小P值,如果小于pEnter,将其纳入selected % 第三步:基于当前selected重新拟合多因素模型,剔除P值大于pRemove的变量 % 第四步:如果本轮没有任何变量进入或剔除,循环终止 changed = stepOnce(X, y, selected, candidateIdx, pEnter, pRemove); if ~changed break; end end % 最后对selected变量拟合一次完整模型,输出系数、OR和置信区间 [coef, stats] = fitFinalModel(X(:, selected), y); end逻辑说明:这里的stepOnce和fitFinalModel是模拟实际代码中的两个内部步骤,压缩包里对应的是一系列MATLAB统计函数调用。循环的终止条件是"某轮迭代中既没有新变量进入、也没有老变量被剔除",这是逐步回归的标准收敛规则。pEnter和pRemove构成一个双向门槛:pEnter管进来的条件,pRemove管出去的条件。用0.05和0.10这两个经典值的意图是让筛选有适度冗余,避免变量在边界上反复进出。
参数说明:如果你的样本量不大(比如事件组只有40-50例),单因素回归的P值会很不稳定,建议把pEnter放宽到0.10、pRemove放到0.15,先让变量进来再说,后面用AIC或BIC去平衡复杂度。另外要注意共线性问题不会直接体现在P值上——如果两个候选变量相关性超过0.7,建议预先只保留其中一个,否则可能出现系数符号异常或标准误爆炸。
3.3 训练集与验证集的切分纪律:origin.xlsx和test.xlsx不能互相偷看
压缩包里已经预置了origin.xlsx(建模数据)和test.xlsx(验证数据),这是一个比随机拆分更干净的数据结构。实际执行时,我的流程是:只在origin.xlsx上做变量筛选和模型拟合,test.xlsx只在最后评估阶段被读入一次。LR_Pz_choice.m、Roc.m的模型参数都只能由训练集决定。有人图省事,先在全量数据上筛变量,再切分训练/测试,这会让验证集的AUC虚高,审稿人一旦追问变量筛选细节就难以解释。
还有一个容易被忽略的点:跑测试集之前,先对比两个数据集的结局事件率。我见过origin事件率30%、test只有8%的情况,这种验证集得到的AUC参考意义很有限,解读时必须在报告里说明事件率差异。正确的做法是检查两个数据集的基线特征表:年龄、性别、关键预测因子的均值与分布是否大致一致,再做结论。这套"验证集先看分布、再算AUC"的动作,能帮你避免很多解释上的尴尬。
跑完LR_Pz_choice.m之后,一定要把最终入选变量的OR值和置信区间打印出来。我一般会在命令行留一份、导出到Excel一份。这个表里的"OR(95%CI)"就是风险预测模型最核心的可交付成果,很多人只盯着AUC而忘了把OR表整理出来,这是本末倒置。
4. ROC曲线与AUC对比:Roc.m、CIAUC.m到DeLong检验的完整链路
4.1 Roc.m:一条ROC曲线从计算到出图
Roc.m是这份资源里被提到最多的脚本,它的输入是预测概率和真实标签,输出是AUC、置信区间以及一张可以直接放进报告里的图。
% Roc.m —— 单模型ROC曲线绘制与AUC计算 % 输入: score为预测概率列向量,y为真实标签(1=事件, 0=非事件) % 输出: auc为曲线下面积,ci为95%置信区间 function [auc, ci] = Roc(score, y) if length(score) ~= length(y) error('score与y长度不一致,请先检查样本对齐'); end % perfcurve是MATLAB统计工具箱的标准ROC函数 [fpr, tpr, ~, auc] = perfcurve(y, score, 1); % AUC标准误的近似计算,事件组较少时改用手动bootstrap n1 = sum(y == 1); n0 = sum(y == 0); seAUC = sqrt(auc * (1 - auc) / min(n1, n0)); ci = [auc - 1.96 * seAUC, auc + 1.96 * seAUC]; figure('Color', 'w'); plot(fpr, tpr, 'LineWidth', 1.8); hold on; plot([0 1], [0 1], 'k--'); xlabel('1 - Specificity'); ylabel('Sensitivity'); title(sprintf('ROC Curve (AUC = %.3f, 95%% CI %.3f - %.3f)', auc, ci(1), ci(2))); legend('Model', 'Random', 'Location', 'SouthEast'); grid on; end逻辑说明:perfcurve会按所有可能的截断阈值计算真阳性率和假阳性率,画出一条从右上角到左下角的阶梯曲线。AUC是所有阈值下判别能力的综合表现:0.5代表随机,0.7-0.8在临床模型里算中等可用,0.9以上要警惕过拟合。代码里用n1和n0中的较小值来近似计算AUC标准误,这个公式来自正态近似,样本量大时够用。
参数说明:score必须是连续预测概率,直接用逻辑回归输出的概率即可;1表示正类(事件)。如果你的事件样本量少于50,这个近似置信区间会偏窄,建议改用bootstrap或者压缩包里的CIAUC.m。ROC test.m则是Roc.m的训练/验证封装版:依次对origin.xlsx和test.xlsx调用Roc.m,把两条曲线画在同一张图上对比,这是论文中常见的呈现方式。
4.2 CIAUC.m与AUC_compare_correlated.m:置信区间与两个模型的差异检验
CIAUC.m干的是"把AUC的置信区间算准"这件事。Roc.m里的近似公式只适用于样本量均衡的情况,CIAUC.m一般会采用DeLong法,它利用排序秩的协方差结构计算AUC方差,比正态近似更稳定,也更容易被审稿人接受。如果你的两组样本量差异很大,直接用Roc.m里的近似区间会有偏差,这时候CIAUC.m就是那个正确答案。
AUC_compare_correlated.m解决的问题更实际:两个模型在同一个测试集上跑出来的AUC分别是0.78和0.82,这个0.04的差异是真实进步还是随机波动?
% AUC_compare_correlated.m —— 两个相关模型AUC差异的DeLong检验 % 输入: y为真实标签,score1、score2为两个模型对同一批样本的预测概率 % 输出: z为检验统计量,p为显著性P值 function [z, p] = AUC_compare_correlated(y, score1, score2) [~, ~, ~, auc1] = perfcurve(y, score1, 1); [~, ~, ~, auc2] = perfcurve(y, score2, 1); % DeLong法核心:计算两个AUC的方差-协方差矩阵 V = delongVariance(y, score1, score2); % 2x2协方差矩阵 d = auc1 - auc2; se = sqrt(V(1,1) + V(2,2) - 2 * V(1,2)); z = d / se; p = 2 * (1 - normcdf(abs(z))); fprintf('AUC1=%.3f, AUC2=%.3f, Delta=%.3f, P=%.4f\n', auc1, auc2, d, p); end逻辑说明:两个模型在相同样本上评估时,它们的AUC并不独立——同一个样本被模型A分错时,很可能也被模型B分错。DeLong法通过计算排序秩的协方差来刻画这种相关性,再检验差值是否显著区别于0,这就是它比朴素z检验稳健的原因。delongVariance对应的是CIAUC.m里按DeLong法构造协方差矩阵的实现,我在这里用函数名做示意。
参数说明:score1和score2必须来自同一批样本,且顺序与y完全一致。我见过有人把训练集的score1和验证集的score2塞进这个函数,算出来的p值极其显著,实际上是不同样本在比,结论完全站不住脚。调用之前,先确认三个变量length一致,再确认score1和score2在同一行上对应的是同一个个体。
4.3 roc.xlsx里到底要导出什么:AUC、最优阈值、灵敏度与特异度
压缩包里的roc.xlsx是Roc.m运行后留下的中间结果,通常会包含每个样本的预测概率、不同阈值下的灵敏度和特异度。写报告或论文时,你需要从中提取三组核心信息:AUC及95%置信区间;Youden指数对应的最优阈值,以及该阈值下的灵敏度与特异度;如果要画校准图,还需要每个十分位组的预测概率均值与实际事件率。
最优阈值一般按Youden指数确定,即灵敏度+特异度-1最大化的点。这个值不一定是0.5,在事件率只有15-20%的数据里,最优阈值经常落在0.2-0.3附近。提取这些指标时别手动读图,写一个循环把结果自动追加到roc.xlsx对应的sheet里,后续做多模型对比会省很多事。另外,这套资源里还能计算PR曲线——当事件率很低时,PR曲线比ROC更能暴露模型在少数类上的表现差异,报告中通常主放ROC、补充PR,这个组合在审稿人那里更耐看。
5. 避坑与常见问题:跑这套预测模型最容易翻车的五个地方
跑风险预测模型最让人头疼的不是统计公式,而是一些看起来不起眼的小问题。以下是压缩包在真实环境下使用时的五个高频翻车现场,按"现象→原因→解决"的格式整理。
5.1 现象:模型AUC高达0.95,验证集上一塌糊涂
原因:大概率是数据泄露。最常见的情况是X矩阵里混入了y那一列本身,或者之前有人用全量数据做过变量筛选再拆分训练/测试,使得验证集的信息已经参与了特征选择。另外,行删除后没有重新对齐标签也会造成类似的假高。
解决:检查LR_Pz_choice.m传入的X是否包含结局变量列;确认所有变量筛选只在origin.xlsx上完成;把test.xlsx单独跑一遍Roc.m,对比训练AUC与验证AUC的差距。如果验证AUC比训练低出0.1以上,优先怀疑数据泄露而不是模型过拟合。快速排查可以用下面这段:
% 快速排查训练矩阵中是否混入结局列 leak = any(all(X == y, 1)); if leak warning('发现疑似结局变量泄漏,请检查变量名'); end上面这段用all(X == y, 1)逐列比较X和y,如果某列对每一行都与y相同,这列极可能就是结局变量本身错放进了预测因子矩阵。这个检查只适用于数值型矩阵,y是逻辑型时先double(y)再比较。
5.2 现象:NRI算出来是负数,但新模型的AUC明明略高
原因:NRI衡量的是"重分类"是否正确,与AUC衡量的整体排序能力是两回事。新模型确实可能把一批人从低风险正确升到高风险,同时也把另一批人错误降级,如果后者净效果更明显,NRI就会为负而AUC反而微升。
解决:先确认NRI.m里的y定义与Roc.m一致(谁是事件组);再看cutoff阈值是否与临床常识匹配。如果事件率约15%却把cutoff设在0.5,绝大多数样本都被分到低风险组,NRI几乎失去区分度。这时候用Category_Free_NRI.m做连续NRI,它会逐一对每个样本的预测概率变化计分,不依赖阈值,结果更稳。
5.3 现象:Roc.m画出来的曲线是锯齿状,不是平滑阶梯
原因:score不是连续的预测概率。要么是逻辑回归输出的概率被四舍五入成了整数评分,要么是score没有按降序排列就直接传给perfcurve。
解决:先检查unique(score)的长度,如果不到10个值,说明概率被离散化,需要检查数据读入时的精度转换;再把score和y按score降序同步排序后传给Roc.m:
% 画ROC前先完成同步排序 [scoreSorted, idx] = sort(score, 'descend'); ySorted = y(idx);默认情况下perfcurve会自己处理排序,但如果你需要手动画阶梯图或者叠加多条曲线,就必须让score和y按同一个idx重排,否则每条曲线对应的样本顺序都不一样,图形就乱了。idx是sort返回的索引,y必须用这个idx重排,这是最关键的一步。
5.4 现象:xlsread读进来后整列是NaN,但Excel里肉眼看着全是数字
原因:这一列的单元格被Excel设成了文本格式,部分单元格左上角还有绿色角标。MATLAB读入时默认把这些当字符处理,转数值时失败就成了NaN。
解决:在Excel里选中整列,用"数据-分列-完成"强制转成常规格式;更省事的办法是让读入脚本用TreatAsMissing选项统一处理缺失标记后再读。清洗完立即检查每列缺失比例:
% 读入后立即检查每一列的缺失比例 nanRatio = mean(isnan(X), 1); disp([varNames; num2cell(nanRatio)]);一行命令就能输出每个变量的NaN占比,如果看到某列是1.0,说明这一列根本没有读入数值。X是数值矩阵,varNames是表头变量名元胞数组,两个对齐后打印出来,比逐个查看Excel省时得多。
5.5 现象:AUC_compare_correlated.m报"矩阵维度不一致",但三个变量看起来都正常
原因:y、score1、score2的行数不一致。最常见的源头是读取Excel时用了不同的range分别读三列,其中一列的行数因为表尾空行被截掉了。
解决:别分开读,直接用readtable把整张表读进来,再用列名取数,这样三个向量天然对齐。如果必须从不同文件取数据,读完立刻用assert(isequal(length(y), length(score1), length(score2)))做一次性校验,比后面报错时再来排查快得多。
补充一条环境问题:在Windows解压这份rar后,__MACOSX和以._开头的文件不会自己消失,MATLAB扫描当前目录时会把它们识别为无扩展名文件,偶尔干扰addpath。解压后先手动清理掉,再打开main.m,能避免很多奇怪的路径报错。
6. NRI快速验证:AUC没涨时,怎么证明新模型仍然有价值
6.1 场景先行:AUC只涨0.01,但不代表模型没有进步
一个常见场景:基线模型AUC=0.78,加入新标志物后AUC=0.79,DeLong检验P=0.20,看起来新变量没用。但如果你逐个样本去看重分类,会发现一批原本被分到中风险的患者被正确提升到了高风险组,这正是AUC这种"整体排序指标"捕捉不到的局部改善。NRI就是为这种场景设计的:它不做全局排序比较,而是直接统计新旧模型对每个样本的风险分类变化,净有多少人从分错变成分对。
6.2 标准NRI与分类自由NRI的选型和使用
% NRI.m —— 标准NRI,需要指定风险分层阈值 % 输入: y为二分类结局,pOld/pNew为旧/新模型对同一批样本的预测概率 % cutoff为风险分层切点向量,例如[0.2, 0.8] [nriStandard, comp] = NRI(y, pOld, pNew, [0.2, 0.8]); % Category_Free_NRI.m —— 分类自由NRI,不需要阈值,适合连续评分比较 [nriFree, compFree] = Category_Free_NRI(y, pOld, pNew);逻辑说明:标准NRI把预测概率切分成低、中、高三段,然后分别统计事件组和非事件组中正确重分类与错误重分类的人数差。分类自由NRI则不做切段,它对每个样本按新旧概率差直接打分,事件组的概率上升记为+1、下降记为-1,非事件组相反,最后求和。两者计算方式不同,但回答的是同一类问题:新模型到底有没有把人群分得更准。
参数说明:pOld和pNew必须来自同一批样本,且顺序与y严格对齐。cutoff的选择要基于临床意义,不要用数据驱动的方式寻优,比如"低风险<20%、高风险>80%"这种定义。还有一个小前置动作:先跑FindNandD.m确认事件组和非事件组人数N和D,如果事件组样本太少,标准NRI的分组会很稀疏,直接改用连续NRI更稳妥。
从那以后,我每次完成一个风险预测模型都会强制走一遍NRI检查:先看AUC和置信区间,再算NRI,两者互相印证。如果AUC变化不显著但NRI显著为正,报告里就重点讲重分类价值;如果NRI为负而AUC上升,我会回头看阈值设置和样本分布。希望这个习惯也能帮你减少一些被审稿人追问的尴尬。
本文还有配套的精品资源,点击获取